]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
updated dphi correlation analysis
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Apr 2011 11:46:08 +0000 (11:46 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Apr 2011 11:46:08 +0000 (11:46 +0000)
added new container class
introduced corr fct calculation as fct of z vtx

15 files changed:
PWG4/CMakelibPWG4JetTasks.pkg
PWG4/JetTasks/AliAnalyseLeadingTrackUE.cxx
PWG4/JetTasks/AliAnalysisTaskLeadingTrackUE.cxx
PWG4/JetTasks/AliAnalysisTaskPhiCorrelations.cxx
PWG4/JetTasks/AliAnalysisTaskPhiCorrelations.h
PWG4/JetTasks/AliTHn.cxx [new file with mode: 0644]
PWG4/JetTasks/AliTHn.h [new file with mode: 0644]
PWG4/JetTasks/AliUEHist.cxx
PWG4/JetTasks/AliUEHist.h
PWG4/JetTasks/AliUEHistograms.cxx
PWG4/JetTasks/AliUEHistograms.h
PWG4/PWG4JetTasksLinkDef.h
PWG4/macros/AddTaskPhiCorrelations.C
PWG4/macros/AddTaskPhiCorrelationsQA.C
PWG4/macros/dphicorrelations/correct.C

index 12c394a4519f7ebaa67d468bb1eae9bb5223e954..af22ab636975878496e3d27763de5726f24e8bc6 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS  JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliHistogramsUE.cxx JetTasks/AliAnalyseUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisTaskJetSpectrum2.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnalysisTaskJetServices.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx JetTasks/AliAnalysisTaskJetCorrections.cxx JetTasks/AliAnalysisTaskThreeJets.cxx JetTasks/AliPWG4HighPtQATPConly.cxx JetTasks/AliPWG4HighPtQAMC.cxx JetTasks/AliPWG4HighPtSpectra.cxx JetTasks/AliPWG4CosmicCandidates.cxx JetTasks/AliAnalysisTaskJetChem.cxx JetTasks/AliAnalysisTaskFragmentationFunction.cxx JetTasks/AliAnalysisTaskMinijet.cxx JetTasks/AliUEHist.cxx JetTasks/AliUEHistograms.cxx JetTasks/AliAnalyseLeadingTrackUE.cxx JetTasks/AliAnalysisTaskLeadingTrackUE.cxx JetTasks/AliAnalysisTaskQGSep.cxx JetTasks/AliAnalysisTaskJetsTM.cxx JetTasks/AliAnalysisTaskJetResponse.cxx JetTasks/AliPWG4HighPtTrackQA.cxx JetTasks/AliEventPoolManager.cxx JetTasks/AliPhiCorrelationsQATask.cxx JetTasks/AliAnalysisTaskPhiCorrelations.cxx JetTasks/AliAnalysisTaskJetResponseV2.cxx)
+set ( SRCS  JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliHistogramsUE.cxx JetTasks/AliAnalyseUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisTaskJetSpectrum2.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnalysisTaskJetServices.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx JetTasks/AliAnalysisTaskJFSystematics.cxx JetTasks/AliAnalysisTaskJetCorrections.cxx JetTasks/AliAnalysisTaskThreeJets.cxx JetTasks/AliPWG4HighPtQATPConly.cxx JetTasks/AliPWG4HighPtQAMC.cxx JetTasks/AliPWG4HighPtSpectra.cxx JetTasks/AliPWG4CosmicCandidates.cxx JetTasks/AliAnalysisTaskJetChem.cxx JetTasks/AliAnalysisTaskFragmentationFunction.cxx JetTasks/AliAnalysisTaskMinijet.cxx JetTasks/AliUEHist.cxx JetTasks/AliUEHistograms.cxx JetTasks/AliAnalyseLeadingTrackUE.cxx JetTasks/AliAnalysisTaskLeadingTrackUE.cxx JetTasks/AliAnalysisTaskQGSep.cxx JetTasks/AliAnalysisTaskJetsTM.cxx JetTasks/AliAnalysisTaskJetResponse.cxx JetTasks/AliPWG4HighPtTrackQA.cxx JetTasks/AliEventPoolManager.cxx JetTasks/AliPhiCorrelationsQATask.cxx JetTasks/AliAnalysisTaskPhiCorrelations.cxx JetTasks/AliAnalysisTaskJetResponseV2.cxx JetTasks/AliTHn.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index ac08e24a02a7b035e9b22afedcdfb7c8276e58e1..635c1d1b37fa9a03087217b1c67fd7f19c3d44aa 100644 (file)
@@ -104,8 +104,8 @@ Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
 {
   
   // select track according to set of cuts
-  if (! fEsdTrackCuts->IsSelected(track) )return kFALSE;
-  if (!fEsdTrackCutsSPD->IsSelected(track) && fEsdTrackCutsSDD->IsSelected(track)) return kFALSE;
+  if (!fEsdTrackCuts->IsSelected(track) )return kFALSE;
+  if (fEsdTrackCutsSPD && fEsdTrackCutsSDD && !fEsdTrackCutsSPD->IsSelected(track) && fEsdTrackCutsSDD->IsSelected(track)) return kFALSE;
 
   return kTRUE;
 }
@@ -114,22 +114,29 @@ Bool_t AliAnalyseLeadingTrackUE::ApplyCuts(TObject* track)
 //____________________________________________________________________
 void AliAnalyseLeadingTrackUE::DefineESDCuts(Int_t /*filterbit*/){
   
-   // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
-   // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
+  // Reproduces the cuts of the corresponding bit in the ESD->AOD filtering
+  // (see $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C)
 
-   fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
-   fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
-   // Add SPD requirement 
-   fEsdTrackCutsSPD = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
-   fEsdTrackCutsSPD->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
-   // Add SDD requirement 
-   fEsdTrackCutsSDD = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
-   fEsdTrackCutsSDD->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
+  if (fFilterBit == 128)
+  {
+    fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+    fEsdTrackCuts->SetMinNClustersTPC(70);
+  }
+  else
+  {
+    fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
+    fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
+
+    // Add SPD requirement 
+    fEsdTrackCutsSPD = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
+    fEsdTrackCutsSPD->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+
+    // Add SDD requirement 
+    fEsdTrackCutsSDD = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");
+    fEsdTrackCutsSDD->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);
+  }
 }
 
-
 //____________________________________________________________________
 TObjArray*  AliAnalyseLeadingTrackUE::FindLeadingObjects(TObject *obj)
 {
@@ -171,6 +178,10 @@ TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject*
 
   Int_t nTracks = NParticles(obj);
   TObjArray* tracks = new TObjArray;
+  
+  // for TPC only tracks
+  if (fFilterBit == 128 && obj->InheritsFrom("AliESDEvent"))
+    tracks->SetOwner(kTRUE);
  
   // Loop over tracks or jets
   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
@@ -379,8 +390,36 @@ AliVParticle*  AliAnalyseLeadingTrackUE::ParticleWithCuts(TObject* obj, Int_t ip
        if (!( ApplyCuts(part)) )
         return 0; 
        
-       // only primary candidates (does not exist for ESD tracks??????)
-       //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
+       if (fFilterBit == 128)
+       {
+         // create TPC only tracks constrained to the SPD vertex
+
+         const AliESDVertex *vtxSPD = esdEvent->GetPrimaryVertexSPD();
+
+         AliESDtrack* track = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent, ipart);
+         if(!track) return 0;
+    
+         // laser warm up tracks
+         if (track->GetTPCsignal() < 10.)
+           return 0;
+    
+         if(track->Pt()>0.){
+           // only constrain tracks above threshold
+           AliExternalTrackParam exParam;
+           // take the B-feild from the ESD, no 3D fieldMap available at this point
+           Bool_t relate = kFALSE;
+           relate = track->RelateToVertex(vtxSPD,esdEvent->GetMagneticField(),kVeryBig,&exParam);
+           if(!relate)
+           {
+//                 Printf("relating failed");
+             delete track;
+             return 0;
+           }
+           track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
+         }
+         
+         part = track;
+       }
        
        // eventually only hadrons
        //TO-DO
@@ -532,12 +571,9 @@ Bool_t  AliAnalyseLeadingTrackUE::TriggerSelection(const TObject* obj)
   if (!obj) // MC
     return kFALSE;
 
-  //Use AliPhysicsSelection to select good events
-  if( !obj->InheritsFrom("AliAODInputHandler") ) { // not needed for AOD input 
-       if (!(((AliInputEventHandler*)obj)->IsEventSelected()&AliVEvent::kMB))return kFALSE;
-        }                                
-
-  // TODO for AOD input trigger bit needs to be checked too (is stored in the AOD)
+  // Use AliPhysicsSelection to select good events, works for ESD and AOD
+  if (!(((AliInputEventHandler*)obj)->IsEventSelected()&(AliVEvent::kMB|AliVEvent::kUserDefined)))
+    return kFALSE;
 
   return kTRUE;
 }
index c09fc9aa2cdc3afca652d5f851f11383710da47f..d55ca87d6d113c13b52efdf3930572f24bcc4da3 100644 (file)
@@ -352,7 +352,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
   // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
   // (MC-true leading particle and MC-true all particles)
   // STEP 0
-  fHistosUE->Fill(fillId,AliUEHist::kCFStepAll,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
+  fHistosUE->Fill(fillId,0,AliUEHist::kCFStepAll,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
   
   // Trigger selection ************************************************
   if (fAnalyseUE->TriggerSelection(fInputHandler))
@@ -370,7 +370,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
            // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
            // (MC-true leading particle and MC-true all particles)
            // STEP 1
-           fHistosUE->Fill(fillId,AliUEHist::kCFStepTriggered,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
+           fHistosUE->Fill(fillId,0,AliUEHist::kCFStepTriggered,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
   
            // count number of MC tracks above 150 MeV/c
            Int_t nMCTracks = 0; 
@@ -391,7 +391,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
              // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
              // (MC-true leading particle and MC-true all particles)
              // STEP 2
-             fHistosUE->Fill(fillId,AliUEHist::kCFStepVertex,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
+             fHistosUE->Fill(fillId,0,AliUEHist::kCFStepVertex,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
     
               // fill here for tracking efficiency
               // loop over particle species
@@ -425,7 +425,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
                // (MC-true leading particle and MC-true all particles)
                // STEP 3
                if (!fReduceMemoryFootprint)
-                 fHistosUE->Fill(fillId,AliUEHist::kCFStepAnaTopology,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
+                 fHistosUE->Fill(fillId,0,AliUEHist::kCFStepAnaTopology,leadingMC,(TList*)regionSortedParticlesMC->At(0),(TList*)regionSortedParticlesMC->At(1),(TList*)regionsMinMaxMC->At(0),(TList*)regionsMinMaxMC->At(1));
       
                //Sort RECO particles w.r.t. MC-leading and return matched (primary) MC particle 
                // (you cannot sort tracks w.r.t. RECO-leading and plot it vs. MC-leading ...)
@@ -435,7 +435,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
                // (MC leading particle and RECO-matched (quantities from MC particle)  all particles)
                // STEP 4
                //if (!fReduceMemoryFootprint)
-                 fHistosUE->Fill(fillId,AliUEHist::kCFStepTrackedOnlyPrim,leadingMC,(TList*)regionSortedParticlesRECOLTMC->At(0),(TList*)regionSortedParticlesRECOLTMC->At(1),(TList*)regionsMinMaxRECOLTMC->At(0),(TList*)regionsMinMaxRECOLTMC->At(1));
+                 fHistosUE->Fill(fillId,0,AliUEHist::kCFStepTrackedOnlyPrim,leadingMC,(TList*)regionSortedParticlesRECOLTMC->At(0),(TList*)regionSortedParticlesRECOLTMC->At(1),(TList*)regionsMinMaxRECOLTMC->At(0),(TList*)regionsMinMaxRECOLTMC->At(1));
                // comparing this step with step 3 (for all-tracks observables) you get the tracking efficiency
         
                //Sort RECO particles w.r.t. MC-leading and return matched (primary+secondary) MC particle 
@@ -447,7 +447,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
                // (MC leading particle and RECO-matched (quantities from MC particle)  all particles)
                // STEP 5
                //if (!fReduceMemoryFootprint)
-                  fHistosUE->Fill(fillId,AliUEHist::kCFStepTracked,leadingMC,(TList*)regionSortedParticlesRECOLTMC2->At(0),(TList*)regionSortedParticlesRECOLTMC2->At(1),(TList*)regionsMinMaxRECOLTMC2->At(0),(TList*)regionsMinMaxRECOLTMC2->At(1));
+                  fHistosUE->Fill(fillId,0,AliUEHist::kCFStepTracked,leadingMC,(TList*)regionSortedParticlesRECOLTMC2->At(0),(TList*)regionSortedParticlesRECOLTMC2->At(1),(TList*)regionsMinMaxRECOLTMC2->At(0),(TList*)regionsMinMaxRECOLTMC2->At(1));
                // comparing this step with step 3 (for all-tracks observables) you get the tracking efficiency
           
                // SWITCH TO RECONSTRUCTED TRACKS  ************************************
@@ -458,7 +458,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
                // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
                // (RECO leading particle and RECO  all particles)
                // STEP 6
-               fHistosUE->Fill(fillId,AliUEHist::kCFStepReconstructed,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
+               fHistosUE->Fill(fillId,0,AliUEHist::kCFStepReconstructed,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
         
                // STEP 8 for reduced efficiency study
                FillReducedEfficiency(fillId, AliUEHist::kCFStepBiasStudy, ltRECO, kFALSE);
@@ -488,7 +488,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseCorrectionMode()
                  // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
                  // (RECO-MATCHED leading particle and RECO all particles)
                  // STEP 7
-                 fHistosUE->Fill(fillId,AliUEHist::kCFStepRealLeading,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
+                 fHistosUE->Fill(fillId,0,AliUEHist::kCFStepRealLeading,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
                  // comparing this step with step 6 (for leading-track observables) you get the efficiency to reconstruct the leading track
                  // comparing this step with step 6 (for all-tracks observables) you see how leading-track misidentification affects the final distributions
                }
@@ -562,7 +562,7 @@ void  AliAnalysisTaskLeadingTrackUE::AnalyseDataMode()
   // Fill UE containers (step, leading track, towards particles, away particles, transverse MIN and MAX particles)
   // (RECO leading particle and RECO  all particles)
   // STEP 6
-  fHistosUE->Fill(eventId,AliUEHist::kCFStepReconstructed,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
+  fHistosUE->Fill(eventId,0,AliUEHist::kCFStepReconstructed,(AliVParticle*)ltRECO->At(0),(TList*)regionSortedParticlesRECO->At(0),(TList*)regionSortedParticlesRECO->At(1),(TList*)regionsMinMaxRECO->At(0),(TList*)regionsMinMaxRECO->At(1));
   
   // STEP 8 and 9 for reduced efficiency study
   FillReducedEfficiency(eventId, AliUEHist::kCFStepBiasStudy, ltRECO, kFALSE);
@@ -619,7 +619,7 @@ void AliAnalysisTaskLeadingTrackUE::FillReducedEfficiency(Int_t eventId, AliUEHi
   TObjArray *regionSortedParticlesRECOLowEff = fAnalyseUE->SortRegions(leading, particleList, 0); 
   TObjArray *regionsMinMaxRECOLowEff = fAnalyseUE->GetMinMaxRegion((TList*)regionSortedParticlesRECOLowEff->At(2), (TList*)regionSortedParticlesRECOLowEff->At(3));
     
-  fHistosUE->Fill(eventId,step,leading,(TList*)regionSortedParticlesRECOLowEff->At(0), (TList*)regionSortedParticlesRECOLowEff->At(1), (TList*)regionsMinMaxRECOLowEff->At(0), (TList*)regionsMinMaxRECOLowEff->At(1));
+  fHistosUE->Fill(eventId,0,step,leading,(TList*)regionSortedParticlesRECOLowEff->At(0), (TList*)regionSortedParticlesRECOLowEff->At(1), (TList*)regionsMinMaxRECOLowEff->At(0), (TList*)regionsMinMaxRECOLowEff->At(1));
 
   delete regionSortedParticlesRECOLowEff;
   delete regionsMinMaxRECOLowEff;
index e7be245bf4adfb3bbfb577b458195dbda5fa5511..3011bd8bfa7f509f04e4b442e70f62c5f6ff36df 100644 (file)
@@ -1,4 +1,3 @@
-
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
@@ -88,7 +87,8 @@ fHistos(0x0),
 fHistosMixed(0),
 fkTrackingEfficiency(0x0),
 // handlers and events
-fAOD(0x0),           
+fAOD(0x0),
+fESD(0x0),
 fArrayMC(0x0),
 fInputHandler(0x0),
 fMcEvent(0x0),
@@ -98,7 +98,7 @@ fPoolMgr(0x0),
 fListOfHistos(0x0), 
 // event QA
 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default) 
-fZVertex(10.),
+fZVertex(7.),
 fCentralityMethod("V0M"),
 // track cuts
 fTrackEtaCut(0.8),
@@ -114,7 +114,6 @@ fSelectCharge(0)
   DefineInput(0, TChain::Class());
   // Output slot #0 writes into a TList container
   DefineOutput(0, TList::Class());
-
 }
 
 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations() 
@@ -140,22 +139,33 @@ void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
   
   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
   
-  if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
-       fAOD = ((AliAODInputHandler*)handler)->GetEvent();
-       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
-  } else {  //output AOD
-       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
-       if( handler && handler->InheritsFrom("AliAODHandler") ) {
-               fAOD = ((AliAODHandler*)handler)->GetAOD();
-               if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
-       } else {  // no AOD
-               AliWarning("I can't get any AOD Event Handler");
-               }
-       }       
+  if( handler && handler->InheritsFrom("AliAODInputHandler") ) 
+  { // input AOD
+    fAOD = ((AliAODInputHandler*)handler)->GetEvent();
+    if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
+  } 
+  else 
+  {  //output AOD
+    handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
+    if (handler && handler->InheritsFrom("AliAODHandler") ) 
+    {
+      fAOD = ((AliAODHandler*)handler)->GetAOD();
+      if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
+    } 
+    else 
+    {  // no AOD
+      AliWarning("I can't get any AOD Event Handler");
+    }
+  }
+  
+  if (handler && handler->InheritsFrom("AliESDInputHandler") ) 
+  { // input ESD
+    // pointer received per event in ::Exec
+    if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
+  } 
   
   // Initialize common pointers
   Initialize();
-   
 }
 
 //____________________________________________________________________
@@ -169,7 +179,7 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
   fAnalyseUE = new AliAnalyseLeadingTrackUE();
   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fPtMin);
   fAnalyseUE->SetDebug(fDebug); 
-  fAnalyseUE->DefineESDCuts(0);
+  fAnalyseUE->DefineESDCuts(fFilterBit);
 
   // Initialize output list of containers
   if (fListOfHistos != NULL){
@@ -200,20 +210,15 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
 
   // event mixing
   //  Int_t trackDepth = 100; // Require e.g. 20 5-track events, or 2 50-track events
-  Int_t trackDepth = 1000; // Require e.g. 20 5-track events, or 2 50-track events
+  Int_t trackDepth = 5000; // Require e.g. 20 5-track events, or 2 50-track events
   Int_t poolsize   = 100;  // Maximum number of events
   
   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
   
-  Int_t nZvtxBins  = 5;
-  Double_t zMin    = -10.0;
-  Double_t zMax    =  10.0;
-  Double_t zVtxBinWidth =  (zMax - zMin) / nZvtxBins;
-  Double_t* zvtxbin = new Double_t[nZvtxBins+1];
-  for (int iz=0; iz<nZvtxBins+1; iz++)
-    zvtxbin[iz] = zMin + iz*zVtxBinWidth;
-  
+  Int_t nZvtxBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
+  Double_t* zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
+
   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
   
   delete[] zvtxbin;
@@ -223,13 +228,24 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
 {
   // array of MC particles
-  if (fMcHandler){
-    fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
-    if (!fArrayMC)
-      AliFatal("No array of MC particles found !!!");
+  if (fMcHandler) {
+    if (fAOD)
+    {
+      fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
+      if (!fArrayMC)
+       AliFatal("No array of MC particles found !!!");
+    }
     fMcEvent = fMcHandler->MCEvent();
   }
 
+  // receive ESD pointer if we are not running AOD analysis
+  if (!fAOD)
+  {
+    AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+    if (handler && handler->InheritsFrom("AliESDInputHandler"))
+      fESD = ((AliESDInputHandler*)handler)->GetEvent();
+  }
+
   // Analyse the event
   if (fMode) AnalyseCorrectionMode();
   else AnalyseDataMode();
@@ -266,7 +282,12 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
   
   if (fCentralityMethod.Length() > 0)
   {
-    AliCentrality *centralityObj = fAOD->GetHeader()->GetCentralityP();
+    AliCentrality *centralityObj = 0;
+    if (fAOD)
+      centralityObj = fAOD->GetHeader()->GetCentralityP();
+    else if (fESD)
+      centralityObj = fESD->GetCentrality();
+    
     if (centralityObj)
     {
       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
@@ -279,19 +300,30 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
      }
   }
   
+  // Support for ESD and AOD based analysis
+  AliVEvent* inputEvent = fAOD;
+  if (!inputEvent)
+    inputEvent = fESD;
+  
+  TObject* mc = fArrayMC;
+  if (!mc)
+    mc = fMcEvent;
+  
   // count all events
   fHistos->FillEvent(centrality, -1);
   
   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
   if (!fAnalyseUE->VertexSelection(fMcEvent, 0, fZVertex)) 
     return;
+  
+  Float_t zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
     
   // Get MC primaries
-  TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(fArrayMC, 0, kTRUE, -1, kTRUE);
+  TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
   
   // (MC-true all particles)
   // STEP 0
-  fHistos->FillCorrelations(centrality, AliUEHist::kCFStepAll, tracksMC);
+  fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC);
   
   // Trigger selection ************************************************
   if (fAnalyseUE->TriggerSelection(fInputHandler))
@@ -299,16 +331,16 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
     fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
     
     // Vertex selection *************************************************
-    if (fAnalyseUE->VertexSelection(fAOD, fnTracksVertex, fZVertex))
+    if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
     {
       // fill here for tracking efficiency
       // loop over particle species
       
       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
       {
-        TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(fArrayMC, 0x0, kTRUE, particleSpecies, kTRUE);
-        TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(fAOD, fArrayMC, kTRUE, particleSpecies, kTRUE);
-        TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(fAOD, fArrayMC, kFALSE, particleSpecies, kTRUE);
+        TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
+        TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
+        TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
       
         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, particleSpecies, centrality);
         
@@ -319,30 +351,31 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
     
       // (MC-true all particles)
       // STEP 2
-      fHistos->FillCorrelations(centrality, AliUEHist::kCFStepVertex, tracksMC);
+      if (!fReduceMemoryFootprint)
+       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC);
       
       // Get MC primaries that match reconstructed track
-      TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(fAOD, fArrayMC, kTRUE, -1, kTRUE);
+      TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
       
       // (RECO-matched (quantities from MC particle) primary particles)
       // STEP 4
-      fHistos->FillCorrelations(centrality, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim);
       
       // Get MC primaries + secondaries that match reconstructed track
-      TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(fAOD, fArrayMC, kFALSE, -1, kTRUE);
+      TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
       
       // (RECO-matched (quantities from MC particle) all particles)
       // STEP 5
-      fHistos->FillCorrelations(centrality, AliUEHist::kCFStepTracked, tracksRecoMatchedAll);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll);
       
       // Get RECO tracks
-      TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(fAOD, 0, kTRUE, -1, kTRUE);
+      TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
       
       // (RECO all tracks)
       // STEP 6
-      fHistos->FillCorrelations(centrality, AliUEHist::kCFStepReconstructed, tracks);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
       
-      if (1)
+      if (0 && !fReduceMemoryFootprint)
       {
         // make list of secondaries (matched with MC)
         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
@@ -351,10 +384,10 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
       
         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
-        fHistos->FillCorrelations(centrality, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll);
+        fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll);
         
         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
-        fHistos->FillCorrelations(centrality, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries);
+        fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries);
       
         // plot delta phi vs process id of secondaries
         // trigger particles: primaries in 4 < pT < 10
@@ -413,15 +446,19 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   if (!fInputHandler)
     return;
     
-  if (!(fInputHandler->IsEventSelected()&AliVEvent::kMB))
+  if (!(fInputHandler->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kUserDefined)))
     return;
 
   Double_t centrality = 0;
   
   if (fCentralityMethod.Length() > 0)
   {
-    //AliCentrality *centralityObj = esd->GetCentrality();
-    AliCentrality *centralityObj = fAOD->GetHeader()->GetCentralityP();
+    AliCentrality *centralityObj = 0;
+    if (fAOD)
+      centralityObj = fAOD->GetHeader()->GetCentralityP();
+    else if (fESD)
+      centralityObj = fESD->GetCentrality();
+    
     if (centralityObj)
       centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
       //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
@@ -430,6 +467,13 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
     Printf("Centrality is %f", centrality);  
   }
   
+  // Support for ESD and AOD based analysis
+  AliVEvent* inputEvent = fAOD;
+  if (!inputEvent)
+    inputEvent = fESD;
+
+  fHistos->SetRunNumber(inputEvent->GetRunNumber());
+  
   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
   
@@ -440,7 +484,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
   
   // Vertex selection *************************************************
-  if(!fAnalyseUE->VertexSelection(fAOD, fnTracksVertex, fZVertex)) return;
+  if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
   
   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
@@ -448,7 +492,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   if (centrality < 0)
     return;
 
-  TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(fAOD, 0, kTRUE, -1, kTRUE);
+  TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
   //Printf("Accepted %d tracks", tracks->GetEntries());
   
   // event mixing
@@ -469,7 +513,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   //    FillCorrelations(). Also nMix should be passed in, so a weight
   //    of 1./nMix can be applied.
 
-  AliAODVertex* vertex = fAOD->GetPrimaryVertex();
+  const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
   Double_t zVtx = vertex->GetZ();
   
   AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
@@ -480,7 +524,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   //pool->SetDebug(1);
     
   // Fill containers at STEP 6 (reconstructed)
-  fHistos->FillCorrelations(centrality, AliUEHist::kCFStepReconstructed, tracks);
+  fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
   
   if (pool->IsReady()) 
   {
@@ -491,7 +535,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
     // Fill mixed-event histos here  
     for (Int_t jMix=0; jMix<nMix; jMix++) {
       TObjArray* bgTracks = pool->GetEvent(jMix);
-      fHistosMixed->FillCorrelations(centrality, AliUEHist::kCFStepReconstructed, tracks, bgTracks, 1.0 / nMix, (jMix == 0));
+      fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, bgTracks, 1.0 / nMix, (jMix == 0));
     }
   }
   
index bd45bbe2bcabf38376fe681f86e823864a1fa5d8..8fa096a2b9cbc94c0b45bf085da66fed3a4627b5 100644 (file)
@@ -43,6 +43,7 @@ class AliVParticle;
 class TH1D;
 class TObjArray;
 class AliEventPoolManager;
+class AliESDEvent;
 
 class  AliAnalysisTaskPhiCorrelations : public AliAnalysisTask
   {
@@ -106,6 +107,7 @@ class  AliAnalysisTaskPhiCorrelations : public AliAnalysisTask
 
     // Handlers and events
     AliAODEvent*             fAOD;             //! AOD Event 
+    AliESDEvent*             fESD;             //! ESD Event 
     TClonesArray*            fArrayMC;         //! Array of MC particles 
     AliInputEventHandler*    fInputHandler;    //! Generic InputEventHandler 
     AliMCEvent*              fMcEvent;         //! MC event
diff --git a/PWG4/JetTasks/AliTHn.cxx b/PWG4/JetTasks/AliTHn.cxx
new file mode 100644 (file)
index 0000000..fb97e97
--- /dev/null
@@ -0,0 +1,346 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, 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.                  *
+ **************************************************************************/
+
+/* $Id: AliTHn.cxx 20164 2007-08-14 15:31:50Z morsch $ */
+
+//
+// this storage container is optimized for small memory usage
+//   under/over flow bins do not exist
+//   sumw2 structure is float only and only create when the weight != 1
+// all histogram functionality (projections, axis, ranges, etc) are taken from THnSparse by propagating 
+// the information up into the THnSparse structure (in the ananalysis *after* data processing and merging)
+//
+// the derivation from THnSparse is obviously against many OO rules. correct would be a common baseclass of THnSparse and THn.
+// 
+// Author: Jan Fiete Grosse-Oetringhaus
+
+#include "AliTHn.h"
+#include "TList.h"
+#include "TCollection.h"
+#include "AliLog.h"
+#include "TArrayF.h"
+#include "THnSparse.h"
+#include "TMath.h"
+
+ClassImp(AliTHn)
+
+AliTHn::AliTHn() : 
+  AliCFContainer(),
+  fNBins(0),
+  fNVars(0),
+  fNSteps(0),
+  fValues(0),
+  fSumw2(0)
+{
+  // Constructor
+}
+
+AliTHn::AliTHn(const Char_t* name, const Char_t* title,const Int_t nSelStep, const Int_t nVarIn, const Int_t* nBinIn) : 
+  AliCFContainer(name, title, nSelStep, nVarIn, nBinIn),
+  fNBins(0),
+  fNVars(nVarIn),
+  fNSteps(nSelStep),
+  fValues(0),
+  fSumw2(0)
+{
+  // Constructor
+
+  fNBins = 1;
+  for (Int_t i=0; i<fNVars; i++)
+    fNBins *= nBinIn[i];
+  
+  Init();
+}
+
+void AliTHn::Init()
+{
+  // initialize
+  
+  fValues = new TArrayF*[fNSteps];
+  fSumw2 = new TArrayF*[fNSteps];
+  
+  for (Int_t i=0; i<fNSteps; i++)
+  {
+    fValues[i] = 0;
+    fSumw2[i] = 0;
+  }
+} 
+
+AliTHn::AliTHn(const AliTHn &c) :
+  AliCFContainer(),
+  fNBins(0),
+  fNVars(0),
+  fNSteps(0),
+  fValues(0),
+  fSumw2(0)
+{
+  //
+  // AliTHn copy constructor
+  //
+
+  ((AliTHn &) c).Copy(*this);
+}
+
+AliTHn::~AliTHn()
+{
+  // Destructor
+  
+  DeleteContainers();
+  
+  if (fValues)
+  {
+    delete fValues;
+    fValues = 0;
+  }
+
+  if (fSumw2)
+  {
+    delete fSumw2;
+    fSumw2 = 0;
+  }
+}
+
+void AliTHn::DeleteContainers()
+{
+  // delete data containers
+  
+  for (Int_t i=0; i<fNSteps; i++)
+  {
+    if (fValues && fValues[i])
+    {
+      delete fValues[i];
+      fValues[i] = 0;
+    }
+    
+    if (fSumw2 && fSumw2[i])
+    {
+      delete fSumw2[i];
+      fSumw2[i] = 0;
+    }
+  }
+}
+
+//____________________________________________________________________
+AliTHn &AliTHn::operator=(const AliTHn &c)
+{
+  // assigment operator
+
+  if (this != &c)
+    ((AliTHn &) c).Copy(*this);
+
+  return *this;
+}
+
+//____________________________________________________________________
+void AliTHn::Copy(TObject& c) const
+{
+  // copy function
+
+  AliTHn& target = (AliTHn &) c;
+  
+  AliCFContainer::Copy(target);
+  
+  target.fNSteps = fNSteps;
+  target.fNBins = fNBins;
+  target.fNVars = fNVars;
+  
+  target.Init();
+
+  for (Int_t i=0; i<fNSteps; i++)
+  {
+    if (fValues[i])
+      target.fValues[i] = new TArrayF(*(fValues[i]));
+    else
+      target.fValues[i] = 0;
+    
+    if (fSumw2[i])
+      target.fSumw2[i] = new TArrayF(*(fSumw2[i]));
+    else
+      target.fSumw2[i] = 0;
+  }
+}
+
+//____________________________________________________________________
+Long64_t AliTHn::Merge(TCollection* list)
+{
+  // Merge a list of AliTHn objects with this (needed for
+  // PROOF). 
+  // Returns the number of merged objects (including this).
+
+  if (!list)
+    return 0;
+  
+  if (list->IsEmpty())
+    return 1;
+  
+  AliCFContainer::Merge(list);
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+  
+  Int_t count = 0;
+  while ((obj = iter->Next())) {
+    
+    AliTHn* entry = dynamic_cast<AliTHn*> (obj);
+    if (entry == 0) 
+      continue;
+
+    for (Int_t i=0; i<fNSteps; i++)
+    {
+      if (entry->fValues[i])
+      {
+       if (!fValues[i])
+         fValues[i] = new TArrayF(fNBins);
+      
+       for (Long64_t l = 0; l<fNBins; l++)
+         fValues[i]->GetArray()[l] += entry->fValues[i]->GetArray()[l];
+      }
+
+      if (entry->fSumw2[i])
+      {
+       if (!fSumw2[i])
+         fSumw2[i] = new TArrayF(fNBins);
+      
+       for (Long64_t l = 0; l<fNBins; l++)
+         fSumw2[i]->GetArray()[l] += entry->fSumw2[i]->GetArray()[l];
+      }
+    }
+    
+    count++;
+  }
+
+  return count+1;
+}
+
+void AliTHn::Fill(const Double_t *var, Int_t istep, Double_t weight)
+{
+  // fills an entry
+
+  // calculate global bin index
+  Long64_t bin = 0;
+  for (Int_t i=0; i<fNVars; i++)
+  {
+    bin *= GetAxis(i, 0)->GetNbins();
+    
+    Int_t tmpBin = GetAxis(i, 0)->FindBin(var[i]);
+//     Printf("%d", tmpBin);
+    // under/overflow not supported
+    if (tmpBin < 1 || tmpBin > GetAxis(i, 0)->GetNbins())
+      return;
+    
+    // bins start from 0 here
+    bin += tmpBin - 1;
+//     Printf("%lld", bin);
+  }
+
+  if (!fValues[istep])
+  {
+    fValues[istep] = new TArrayF(fNBins);
+    AliInfo(Form("Created values container for step %d", istep));
+  }
+
+  if (weight != 1)
+  {
+    // initialize with already filled entries (which have been filled with weight == 1), in this case fSumw2 := fValues
+    if (!fSumw2[istep])
+    {
+      fSumw2[istep] = new TArrayF(*fValues[istep]);
+      AliInfo(Form("Created sumw2 container for step %d", istep));
+    }
+  }
+
+  fValues[istep]->GetArray()[bin] += weight;
+  if (fSumw2[istep])
+    fSumw2[istep]->GetArray()[bin] += weight * weight;
+  
+//   Printf("%f", fValues[istep][bin]);
+  
+  // debug
+//   AliCFContainer::Fill(var, istep, weight);
+}
+
+Long64_t AliTHn::GetGlobalBinIndex(const Int_t* binIdx)
+{
+  // calculates global bin index
+  // binIdx contains TAxis bin indexes
+  // here bin count starts at 0 because we do not have over/underflow bins
+  
+  Long64_t bin = 0;
+  for (Int_t i=0; i<fNVars; i++)
+  {
+    bin *= GetAxis(i, 0)->GetNbins();
+    bin += binIdx[i] - 1;
+  }
+
+  return bin;
+}
+
+void AliTHn::FillParent()
+{
+  // fills the information stored in the buffer in this class into the baseclass containers
+  
+  for (Int_t i=0; i<fNSteps; i++)
+  {
+    if (!fValues[i])
+      continue;
+      
+    Float_t* source = fValues[i]->GetArray();
+    // if fSumw2 is not stored, the sqrt of the number of bin entries in source is filled below; otherwise we use fSumw2
+    Float_t* sourceSumw2 = source;
+    if (fSumw2[i])
+      sourceSumw2 = fSumw2[i]->GetArray();
+    
+    THnSparse* target = GetGrid(i)->GetGrid();
+    
+    Int_t* binIdx = new Int_t[fNVars];
+    for (Int_t j=0; j<fNVars; j++)
+      binIdx[j] = 1;
+    
+    Long64_t count = 0;
+    
+    while (1)
+    {
+//       for (Int_t j=0; j<fNVars; j++)
+//     printf("%d ", binIdx[j]);
+      
+      Long64_t globalBin = GetGlobalBinIndex(binIdx);
+//       Printf(" --> %lld", globalBin);
+      
+      if (source[globalBin] != 0)
+      {
+       target->SetBinContent(binIdx, source[globalBin]);
+       target->SetBinError(binIdx, TMath::Sqrt(sourceSumw2[globalBin]));
+       
+       count++;
+      }
+      
+      binIdx[fNVars-1]++;
+      
+      for (Int_t j=fNVars-1; j>0; j--)
+      {
+       if (binIdx[j] > target->GetAxis(j)->GetNbins())
+       {
+         binIdx[j] = 1;
+         binIdx[j-1]++;
+       }
+      }
+      
+      if (binIdx[0] > target->GetAxis(0)->GetNbins())
+       break;
+    }
+    
+    AliInfo(Form("Step %d: copied %lld entries out of %lld bins", i, count, GetGlobalBinIndex(binIdx)));
+  }
+}
diff --git a/PWG4/JetTasks/AliTHn.h b/PWG4/JetTasks/AliTHn.h
new file mode 100644 (file)
index 0000000..39ecba2
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef AliTHn_H
+#define AliTHn_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// optimized data container reusing functionality of AliCFContainer / THnsparse
+
+#include "TObject.h"
+#include "TString.h"
+#include "AliCFContainer.h"
+
+class TArrayF;
+class TCollection;
+
+class AliTHn : public AliCFContainer
+{
+ public:
+  AliTHn();
+  AliTHn(const Char_t* name, const Char_t* title,const Int_t nSelStep, const Int_t nVarIn, const Int_t* nBinIn);
+  
+  virtual ~AliTHn();
+  
+  virtual void  Fill(const Double_t *var, Int_t istep, Double_t weight=1.) ;
+  virtual void  FillParent();
+  
+  TArrayF* GetValues(Int_t step) { return fValues[step]; }
+  TArrayF* GetSumw2(Int_t step)  { return fSumw2[step]; }
+  
+  void DeleteContainers();
+  
+  AliTHn(const AliTHn &c);
+  AliTHn& operator=(const AliTHn& corr);
+  virtual void Copy(TObject& c) const;
+
+  virtual Long64_t Merge(TCollection* list);
+  
+protected:
+  void Init();
+  Long64_t GetGlobalBinIndex(const Int_t* binIdx);
+  
+  Long64_t fNBins;   // number of total bins
+  Int_t    fNVars;   // number of variables
+  Int_t    fNSteps;  // number of selection steps
+  TArrayF **fValues; //[fNSteps] data container
+  TArrayF **fSumw2;  //[fNSteps] data container
+  
+  ClassDef(AliTHn, 2) // THn like container
+};
+
+#endif
index f7cb29a190a2c606f2b63c2dcb577e3b8a112489..f74935d8b84ed93311e3250592fb2b5b3f91ae90 100644 (file)
@@ -34,6 +34,7 @@
 #include "AliLog.h"
 #include "TCanvas.h"
 #include "TF1.h"
+#include "AliTHn.h"
 
 ClassImp(AliUEHist)
 
@@ -50,6 +51,8 @@ AliUEHist::AliUEHist(const char* reqHist) :
   fPtMax(0),
   fCentralityMin(0),
   fCentralityMax(0),
+  fZVtxMin(0),
+  fZVtxMax(0),
   fContaminationEnhancement(0),
   fCombineMinMax(0),
   fCache(0),
@@ -70,9 +73,9 @@ AliUEHist::AliUEHist(const char* reqHist) :
     
   // track level
   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
-  Int_t iTrackBin[5];
-  Double_t* trackBins[5];
-  const char* trackAxisTitle[5];
+  Int_t iTrackBin[6];
+  Double_t* trackBins[6];
+  const char* trackAxisTitle[6];
   
   // eta
   iTrackBin[0] = 20;
@@ -83,14 +86,14 @@ AliUEHist::AliUEHist(const char* reqHist) :
   trackAxisTitle[0] = "#eta";
   
   // delta eta
-  const Int_t kNDeltaEtaBins = 32;
+  const Int_t kNDeltaEtaBins = 40;
   Double_t deltaEtaBins[kNDeltaEtaBins+1];
   for (Int_t i=0; i<=kNDeltaEtaBins; i++)
-    deltaEtaBins[i] = -1.6 + 0.1 * i;
+    deltaEtaBins[i] = -2.0 + 0.1 * i;
   
   // pT
-  iTrackBin[1] = 39;
-  Double_t pTBins[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
+  iTrackBin[1] = 22;
+  Double_t pTBins[] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
   trackBins[1] = pTBins;
   trackAxisTitle[1] = "p_{T} (GeV/c)";
   
@@ -101,11 +104,11 @@ AliUEHist::AliUEHist(const char* reqHist) :
     leadingpTBins[i] = 0.5 * i;
   
   // pT,lead binning 2
-  const Int_t kNLeadingpTBins2 = 14;
-  Double_t leadingpTBins2[] = { 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 100.0 };
+  const Int_t kNLeadingpTBins2 = 9;
+  Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
   
   // phi,lead; this binning starts at -pi/2 and is modulo 3
-  const Int_t kNLeadingPhiBins = 36;
+  const Int_t kNLeadingPhiBins = 72;
   Double_t leadingPhiBins[kNLeadingPhiBins+1];
   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
     leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiBins * i * TMath::TwoPi();
@@ -117,17 +120,21 @@ AliUEHist::AliUEHist(const char* reqHist) :
     multiplicityBins[i] = -0.5 + i;
   multiplicityBins[kNMultiplicityBins] = 200;
   
+  trackBins[3] = multiplicityBins;
+  iTrackBin[3] = kNMultiplicityBins;
+  trackAxisTitle[3] = "multiplicity";
+  
   // centrality (in %)
-  const Int_t kNCentralityBins = 5 + 3 + 8;
-  Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
+  const Int_t kNCentralityBins = 5 + 1 + 9;
+  Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
   
   // particle species
   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
   
-  trackBins[3] = multiplicityBins;
-  iTrackBin[3] = kNMultiplicityBins;
-  trackAxisTitle[3] = "multiplicity";
+  // vtx-z axis
+  const Int_t kNVertexBins = 7;
+  Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
   
   // selection depending on requested histogram
   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
@@ -178,7 +185,7 @@ AliUEHist::AliUEHist(const char* reqHist) :
   }
   else if (axis == 2)
   {
-    nTrackVars = 5;
+    nTrackVars = 6;
     initRegions = 1;
   
     iTrackBin[0] = kNDeltaEtaBins;
@@ -196,11 +203,18 @@ AliUEHist::AliUEHist(const char* reqHist) :
     iTrackBin[4] = kNLeadingPhiBins;
     trackBins[4] = leadingPhiBins;
     trackAxisTitle[4] = "#Delta#phi (rad.)";
+
+    iTrackBin[5] = kNVertexBins;
+    trackBins[5] = vertexBins;
+    trackAxisTitle[5] = "z-vtx (cm)";
   }
     
   for (UInt_t i=0; i<initRegions; i++)
   {
-    fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
+    if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
+      fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
+    else
+      fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
     
     for (Int_t j=0; j<nTrackVars; j++)
     {
@@ -211,8 +225,22 @@ AliUEHist::AliUEHist(const char* reqHist) :
     SetStepNames(fTrackHist[i]);
   }
   
+  // event level
+  Int_t nEventVars = 2;
+  Int_t iEventBin[3];
+
   // track 3rd and 4th axis --> event 1st and 2nd axis
-  fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
+  iEventBin[0] = iTrackBin[2];
+  iEventBin[1] = iTrackBin[3];
+  
+  // plus track 5th axis (in certain cases)
+  if (axis == 2)
+  {
+    nEventVars = 3;
+    iEventBin[2] = iTrackBin[5];
+  }
+  
+  fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
   
   fEventHist->SetBinLimits(0, trackBins[2]);
   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
@@ -220,6 +248,12 @@ AliUEHist::AliUEHist(const char* reqHist) :
   fEventHist->SetBinLimits(1, trackBins[3]);
   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
   
+  if (axis == 2)
+  {
+    fEventHist->SetBinLimits(2, trackBins[5]);
+    fEventHist->SetVarTitle(2, trackAxisTitle[5]);
+  }
+
   SetStepNames(fEventHist);
   
   iTrackBin[2] = kNSpeciesBins;
@@ -247,6 +281,8 @@ AliUEHist::AliUEHist(const AliUEHist &c) :
   fPtMax(0),
   fCentralityMin(0),
   fCentralityMax(0),
+  fZVtxMin(0),
+  fZVtxMax(0),
   fContaminationEnhancement(0),
   fCombineMinMax(0),
   fCache(0),
@@ -335,6 +371,8 @@ void AliUEHist::Copy(TObject& c) const
   target.fPtMax = fPtMax;
   target.fCentralityMin = fCentralityMin;
   target.fCentralityMax = fCentralityMax;
+  target.fZVtxMin = fZVtxMin;
+  target.fZVtxMax = fZVtxMax;
   
   if (fContaminationEnhancement)
     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
@@ -496,7 +534,7 @@ void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_
 }  
 
 //____________________________________________________________________
-TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD)
+TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm)
 {
   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
   //
@@ -507,12 +545,22 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
   //       1: 2D histogram, deltaphi, deltaeta
   //       10: 1D histogram, within |deltaeta| < 0.8
   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
+  //
+  // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
   
   // unzoom all axes
   ResetBinLimits(fTrackHist[region]->GetGrid(step));
   ResetBinLimits(fEventHist->GetGrid(step));
   
   SetBinLimits(fTrackHist[region]->GetGrid(step));
+  
+  // z vtx
+  if (fZVtxMax > fZVtxMin)
+  {
+    Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
+    fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
+    fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
+  }
     
   TH1D* tracks = 0;
   
@@ -536,7 +584,9 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
     Float_t phiRegion = TMath::TwoPi() / 3;
     if (!fCombineMinMax && region == kMin)
       phiRegion /= 2;
-    Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
+    Float_t normalization = phiRegion;
+    if (etaNorm)
+      normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
     //Printf("Normalization: %f", normalization);
     tracks->Scale(1.0 / normalization);
     
@@ -595,16 +645,20 @@ TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Floa
     // normalize to get a density (deta dphi)
     // TODO normalization may be off for 2d histogram
     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
-    TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
-    if (strcmp(axis->GetTitle(), "#eta") == 0)
-    {
-      Printf("Normalizing using eta-axis range");
-      normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
-    }
-    else
+    
+    if (etaNorm)
     {
-      Printf("Normalizing assuming accepted range of |eta| < 0.8");
-      normalization *= 0.8 * 2;
+      TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
+      if (strcmp(axis->GetTitle(), "#eta") == 0)
+      {
+       Printf("Normalizing using eta-axis range");
+       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
+      }
+      else
+      {
+       Printf("Normalizing assuming accepted range of |eta| < 0.8");
+       normalization *= 0.8 * 2;
+      }
     }
     
     //Printf("Normalization: %f", normalization);
index b53b39a931542a791c9e7425bae30130d15a854f..fdfe16d4c8a4f1ebb3b90f530d54852e8a2cfd4b 100644 (file)
@@ -46,7 +46,7 @@ class AliUEHist : public TObject
   
   void CopyReconstructedData(AliUEHist* from);
   
-  TH1* GetUEHist(CFStep step, Region region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1, Int_t multBinBegin = 0, Int_t multBinEnd = -1, Int_t twoD = 0);
+  TH1* GetUEHist(CFStep step, Region region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1, Int_t multBinBegin = 0, Int_t multBinEnd = -1, Int_t twoD = 0, Bool_t etaNorm = kTRUE);
   TH1* GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization = kFALSE);
   
   TH1* GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2 = -1, Int_t source = 1, Int_t axis3 = -1);
@@ -83,6 +83,7 @@ class AliUEHist : public TObject
   void SetEtaRange(Float_t etaMin, Float_t etaMax) { fEtaMin = etaMin; fEtaMax = etaMax; }
   void SetPtRange(Float_t ptMin, Float_t ptMax)    { fPtMin = ptMin; fPtMax = ptMax; }
   void SetCentralityRange(Float_t min, Float_t max)    { fCentralityMin = min; fCentralityMax = max; }
+  void SetZVtxRange(Float_t min, Float_t max)          { fZVtxMin = min; fZVtxMax = max; }
   
   void SetContaminationEnhancement(TH1F* hist)    { fContaminationEnhancement = hist; }
   
@@ -117,6 +118,8 @@ protected:
   Float_t fPtMax;                     // pT max for projections (for track pT, not pT,lead)
   Float_t fCentralityMin;             // centrality min for projections
   Float_t fCentralityMax;             // centrality max for projections
+  Float_t fZVtxMin;                   // z vtx min for projections
+  Float_t fZVtxMax;                   // z vtx max for projections
   
   TH1F* fContaminationEnhancement;    // histogram that contains the underestimation of secondaries in the MC as function of pT
   
@@ -126,7 +129,7 @@ protected:
   
   TString fHistogramType;             // what is stored in this histogram
   
-  ClassDef(AliUEHist, 7) // underlying event histogram container
+  ClassDef(AliUEHist, 8) // underlying event histogram container
 };
 
 #endif
index f7810392260533dc877df0809430bfbf57015d6c..ca318492343924ca8e8396a99dea92204d6a45dd 100644 (file)
@@ -54,7 +54,8 @@ AliUEHistograms::AliUEHistograms(const char* name, const char* histograms) :
   fVertexContributors(0),
   fCentralityDistribution(0),
   fITSClusterMap(0),
-  fSelectCharge(0)
+  fSelectCharge(0),
+  fRunNumber(0)
 {
   // Constructor
   //
@@ -141,7 +142,8 @@ AliUEHistograms::AliUEHistograms(const AliUEHistograms &c) :
   fVertexContributors(0),
   fCentralityDistribution(0),
   fITSClusterMap(0),
-  fSelectCharge(0)
+  fSelectCharge(0),
+  fRunNumber(0)
 {
   //
   // AliUEHistograms copy constructor
@@ -269,7 +271,7 @@ Int_t AliUEHistograms::CountParticles(TList* list, Float_t ptMin)
 }
   
 //____________________________________________________________________
-void AliUEHistograms::Fill(Int_t eventType, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max)
+void AliUEHistograms::Fill(Int_t eventType, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max)
 {
   // fills the UE event histograms
   //
@@ -290,14 +292,15 @@ void AliUEHistograms::Fill(Int_t eventType, AliUEHist::CFStep step, AliVParticle
     multiplicity += CountParticles(min, ptMin);
     multiplicity += CountParticles(max, ptMin);
      
-    FillRegion(AliUEHist::kToward, step, leading, toward, multiplicity);
-    FillRegion(AliUEHist::kAway,   step, leading, away, multiplicity);
-    FillRegion(AliUEHist::kMin,    step, leading, min, multiplicity);
-    FillRegion(AliUEHist::kMax,    step, leading, max, multiplicity);
+    FillRegion(AliUEHist::kToward, zVtx, step, leading, toward, multiplicity);
+    FillRegion(AliUEHist::kAway,   zVtx, step, leading, away, multiplicity);
+    FillRegion(AliUEHist::kMin,    zVtx, step, leading, min, multiplicity);
+    FillRegion(AliUEHist::kMax,    zVtx, step, leading, max, multiplicity);
  
-    Double_t vars[2];
+    Double_t vars[3];
     vars[0] = leading->Pt();
     vars[1] = multiplicity;
+    vars[2] = zVtx;
     for (Int_t i=0; i<fgkUEHists; i++)
       if (GetUEHist(i))
         GetUEHist(i)->GetEventHist()->Fill(vars, step);
@@ -309,7 +312,7 @@ void AliUEHistograms::Fill(Int_t eventType, AliUEHist::CFStep step, AliVParticle
 }
   
 //____________________________________________________________________
-void AliUEHistograms::FillRegion(AliUEHist::Region region, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity)
+void AliUEHistograms::FillRegion(AliUEHist::Region region, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity)
 {
   // loops over AliVParticles in list and fills the given region at the given step
   //
@@ -319,7 +322,7 @@ void AliUEHistograms::FillRegion(AliUEHist::Region region, AliUEHist::CFStep ste
   {
     AliVParticle* particle = (AliVParticle*) list->At(i);
     
-    Double_t vars[5];
+    Double_t vars[6];
     vars[0] = particle->Eta();
     vars[1] = particle->Pt();
     vars[2] = leading->Pt();
@@ -329,7 +332,8 @@ void AliUEHistograms::FillRegion(AliUEHist::Region region, AliUEHist::CFStep ste
       vars[4] -= TMath::TwoPi();
     if (vars[4] < -0.5 * TMath::Pi())
       vars[4] += TMath::TwoPi();
-
+    vars[5] = zVtx;
+    
     if (fNumberDensitypT)
       fNumberDensitypT->GetTrackHist(region)->Fill(vars, step);
       
@@ -379,7 +383,7 @@ void AliUEHistograms::Fill(AliVParticle* leadingMC, AliVParticle* leadingReco)
 }
 
 //____________________________________________________________________
-void AliUEHistograms::FillCorrelations(Double_t centrality, AliUEHist::CFStep step, TSeqCollection* particles, TSeqCollection* mixed, Float_t weight, Bool_t firstTime)
+void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TSeqCollection* particles, TSeqCollection* mixed, Float_t weight, Bool_t firstTime)
 {
   // fills the fNumberDensityPhi histogram
   //
@@ -403,8 +407,8 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, AliUEHist::CFStep st
         fCorrelationpT->Fill(centrality, triggerParticle->Pt());
         fCorrelationEta->Fill(centrality, triggerParticle->Eta());
         fCorrelationPhi->Fill(centrality, triggerParticle->Phi());
-        if (dynamic_cast<AliAODTrack*>(triggerParticle))
-          fITSClusterMap->Fill(((AliAODTrack*) triggerParticle)->GetITSClusterMap(), centrality, triggerParticle->Pt());
+/*        if (dynamic_cast<AliAODTrack*>(triggerParticle))
+          fITSClusterMap->Fill(((AliAODTrack*) triggerParticle)->GetITSClusterMap(), centrality, triggerParticle->Pt());*/
       }
         
       for (Int_t j=0; j<jMax; j++)
@@ -436,7 +440,7 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, AliUEHist::CFStep st
             continue;
         }
         
-        Double_t vars[5];
+        Double_t vars[6];
         vars[0] = triggerParticle->Eta() - particle->Eta();
         vars[1] = particle->Pt();
         vars[2] = triggerParticle->Pt();
@@ -446,6 +450,7 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, AliUEHist::CFStep st
           vars[4] -= TMath::TwoPi();
         if (vars[4] < -0.5 * TMath::Pi())
           vars[4] += TMath::TwoPi();
+       vars[5] = zVtx;
     
         // fill all in toward region and do not use the other regions
         fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->Fill(vars, step, weight);
@@ -454,9 +459,10 @@ void AliUEHistograms::FillCorrelations(Double_t centrality, AliUEHist::CFStep st
       if (firstTime)
       {
         // once per trigger particle
-        Double_t vars[2];
+        Double_t vars[3];
         vars[0] = triggerParticle->Pt();
         vars[1] = centrality;
+       vars[2] = zVtx;
         fNumberDensityPhi->GetEventHist()->Fill(vars, step);
       }
     }
@@ -540,6 +546,16 @@ void AliUEHistograms::SetPtRange(Float_t ptMin, Float_t ptMax)
       GetUEHist(i)->SetPtRange(ptMin, ptMax);
 }
 
+//____________________________________________________________________
+void AliUEHistograms::SetZVtxRange(Float_t min, Float_t max)
+{
+  // sets pT min and max for all contained AliUEHist classes
+  
+  for (Int_t i=0; i<fgkUEHists; i++)
+    if (GetUEHist(i))
+      GetUEHist(i)->SetZVtxRange(min, max);
+}
+
 //____________________________________________________________________
 void AliUEHistograms::SetContaminationEnhancement(TH1F* hist)
 {
@@ -631,6 +647,7 @@ void AliUEHistograms::Copy(TObject& c) const
     target.fITSClusterMap = dynamic_cast<TH3F*> (fITSClusterMap->Clone());
     
   target.fSelectCharge = fSelectCharge;
+  target.fRunNumber = fRunNumber;
 }
 
 //____________________________________________________________________
index 4115fc996d585783173a68ee945ead81425d1b74..a40004da8fe28852009f96ab9ed085a12cb8f80c 100644 (file)
@@ -25,8 +25,8 @@ class AliUEHistograms : public TNamed
   AliUEHistograms(const char* name = "AliUEHistograms", const char* histograms = "123");
   virtual ~AliUEHistograms();
   
-  void Fill(Int_t eventType, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max);
-  void FillCorrelations(Double_t centrality, AliUEHist::CFStep step, TSeqCollection* particles, TSeqCollection* mixed = 0, Float_t weight = 1, Bool_t firstTime = kTRUE);
+  void Fill(Int_t eventType, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max);
+  void FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TSeqCollection* particles, TSeqCollection* mixed = 0, Float_t weight = 1, Bool_t firstTime = kTRUE);
   void Fill(AliVParticle* leadingMC, AliVParticle* leadingReco);
   void FillEvent(Int_t eventType, Int_t step);
   void FillEvent(Double_t centrality, Int_t step);
@@ -44,6 +44,8 @@ class AliUEHistograms : public TNamed
   void SetSumpT(AliUEHist* obj) { fSumpT = obj; }
   void SetNumberDensityPhi(AliUEHist* obj) { fNumberDensityPhi = obj; }
   
+  void SetRunNumber(Long64_t runNumber) { fRunNumber = runNumber; }
+  
   TH2F* GetCorrelationpT()  { return fCorrelationpT; }
   TH2F* GetCorrelationEta() { return fCorrelationEta; }
   TH2F* GetCorrelationPhi() { return fCorrelationPhi; }
@@ -55,11 +57,13 @@ class AliUEHistograms : public TNamed
   TH3F* GetEventCountDifferential() { return fEventCountDifferential; }
   TH1F* GetVertexContributors() { return fVertexContributors; }
   TH1F* GetCentralityDistribution() { return fCentralityDistribution; }
+  Long64_t GetRunNumber() { return fRunNumber; }
   
   void Correct(AliUEHistograms* corrections);
   
   void SetEtaRange(Float_t etaMin, Float_t etaMax);
   void SetPtRange(Float_t ptMin, Float_t ptMax);
+  void SetZVtxRange(Float_t min, Float_t max);
   void SetContaminationEnhancement(TH1F* hist);
   void SetCombineMinMax(Bool_t flag);
   void SetSelectCharge(Int_t selectCharge) { fSelectCharge = selectCharge; }
@@ -75,7 +79,7 @@ class AliUEHistograms : public TNamed
   void Scale(Double_t factor);
   
 protected:
-  void FillRegion(AliUEHist::Region region, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity);
+  void FillRegion(AliUEHist::Region region, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity);
   Int_t CountParticles(TList* list, Float_t ptMin);
   
   static const Int_t fgkUEHists; // number of histograms
@@ -101,7 +105,9 @@ protected:
   
   Int_t fSelectCharge;           // (un)like sign selection when building correlations: 0: no selection; 1: unlike sign; 2: like sign
   
-  ClassDef(AliUEHistograms, 6)  // underlying event histogram container
+  Long64_t fRunNumber;           // run number that has been processed
+  
+  ClassDef(AliUEHistograms, 7)  // underlying event histogram container
 };
 
 #endif
index 934dd3d251cf1345d7836c6b41261527974484c2..ee38c371245212f80ac6ffff3c6576808a44d1b4 100644 (file)
@@ -42,5 +42,6 @@
 #pragma link C++ class AliEventPool+;
 #pragma link C++ class AliEventPoolManager+;
 #pragma link C++ class AliPhiCorrelationsQATask+;
+#pragma link C++ class AliTHn+;
 
 #endif
index e3b5c8ffe9a4c9611d6af267e9eeddd938f0b1e5..a059590a9ead96dff827dff67c59f98db2283127 100644 (file)
@@ -1,4 +1,4 @@
-AliAnalysisTaskPhiCorrelations *AddTaskPhiCorrelations(Int_t analysisMode = 0)\r
+AliAnalysisTaskPhiCorrelations *AddTaskPhiCorrelations(Int_t analysisMode = 0, Bool_t ppRun = kFALSE)\r
 {\r
   // Get the pointer to the existing analysis manager via the static access method.\r
   //==============================================================================\r
@@ -20,25 +20,33 @@ AliAnalysisTaskPhiCorrelations *AddTaskPhiCorrelations(Int_t analysisMode = 0)
   ana->SetDebugLevel(0); \r
   //  ana->SetFilterBit(16);  \r
   //ana->SetFilterBit(64+32);  \r
-  ana->SetFilterBit(1);  \r
-  ana->SetTrackEtaCut(0.8);\r
+  \r
+//   Int_t bit = 1;\r
+  Int_t bit = 128;\r
+  ana->SetFilterBit(bit);  \r
+  \r
+  Printf("AddTaskPhiCorrelations:\n\n\n++++++++++ Using bit %d ++++++++++++\n\n\n", bit);\r
+  \r
+  ana->SetTrackEtaCut(1.0);\r
   ana->SetPtMin(0.15);\r
   //ana->SetEventSelectionBit(AliAnalysisHelperJetTasks::kIsPileUp);\r
   ana->SetReduceMemoryFootprint(kTRUE);\r
   //ana->SetSelectCharge(2);\r
   \r
-  if (1)\r
+  if (0)\r
   {\r
-    Printf("\n\n\n+++++++++++++++ Configuring for p+p! +++++++++++++++++\n\n\n");\r
-    ana->SetCentralityMethod(""); // only for pp\r
+    Printf("AddTaskPhiCorrelations:\n\n\n++++++++++ Using SPD centrality selection ++++++++++++\n\n\n");\r
+    ana->SetCentralityMethod("CL1");\r
   }    \r
   \r
-  if (0)\r
+  if (ppRun)\r
   {\r
-    Printf("\n\n\n++++++++++ Using SPD centrality selection ++++++++++++\n\n\n");\r
-    ana->SetCentralityMethod("CL1");\r
+    Printf("AddTaskPhiCorrelations:\n\n\n+++++++++++++++ Configuring for p+p! +++++++++++++++++\n\n\n");\r
+    ana->SetCentralityMethod(""); // only for pp\r
   }    \r
   \r
+//   gSystem->Sleep(3000);\r
+  \r
   mgr->AddTask(ana);\r
   \r
   // Create ONLY the output containers for the data produced by the task.\r
index 4dda6960fde43a6cfae12026ffac0e80127d4507..60cfe663bb2b3c98eb04c96e6bdc0c62eb6e4f4c 100644 (file)
@@ -12,7 +12,7 @@ AliPhiCorrelationsQATask *AddTaskPhiCorrelationsQA()
   //===========================================================================\r
   AliPhiCorrelationsQATask* ana = new  AliPhiCorrelationsQATask("");\r
   \r
-  ana->SelectCollisionCandidates();\r
+  ana->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kUserDefined);\r
   \r
   mgr->AddTask(ana);\r
   \r
index 97a3a8c5e3a2130c72b0cf22def76c521508e45f..eba2355ab90947392dc412d07e83f39baffbc95e 100644 (file)
@@ -324,29 +324,49 @@ void loadlibs()
   gSystem->Load("libPWG4JetTasks");
 }
 
+const char* lastFileName = 0;
+void* cacheSameEvent = 0;
+void* cacheMixedEvent = 0;
+
 void* GetUEHistogram(const char* fileName, TList** listRef = 0, Bool_t mixed = kFALSE)
 {
-  file = TFile::Open(fileName);
-  if (!file)
-    return 0;
-  
-  list = (TList*) gFile->Get("PWG4_LeadingTrackUE/histosLeadingTrackUE");
-  if (!list)
-    list = (TList*) gFile->Get("PWG4_PhiCorrelations/histosPhiCorrelations");
-    
-  if (!list)
-    return 0;
+  if (!lastFileName || strcmp(lastFileName, fileName) != 0)
+  {
+    lastFileName = fileName;
+    file = TFile::Open(fileName);
+    if (!file)
+      return 0;
     
-  if (listRef)
-    *listRef = list;
+    list = (TList*) gFile->Get("PWG4_LeadingTrackUE/histosLeadingTrackUE");
+    if (!list)
+      list = (TList*) gFile->Get("PWG4_PhiCorrelations/histosPhiCorrelations");
+      
+    if (!list)
+      return 0;
+      
+    if (listRef)
+      *listRef = list;
+      
+    cacheMixedEvent = list->FindObject("AliUEHistogramsMixed");
+    cacheSameEvent = list->FindObject("AliUEHistogramsSame");
+
+    if (mixed)
+      return cacheMixedEvent;
     
-  if (mixed)
-    return list->FindObject("AliUEHistogramsMixed");
-  
-  if (list->FindObject("AliUEHistograms"))
-    return list->FindObject("AliUEHistograms");
+    if (list->FindObject("AliUEHistograms"))
+      return list->FindObject("AliUEHistograms");
+      
+    return cacheSameEvent;
+  }
+  else
+  {
+    Printf("GetUEHistogram --> Using cache");
     
-  return list->FindObject("AliUEHistogramsSame");
+    if (mixed)
+      return cacheMixedEvent;
+    else
+      return cacheSameEvent;
+  }
 }
 
 void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region, Int_t ueHist)
@@ -1275,6 +1295,196 @@ TGraph* GetFlow05()
   return graph;
 }
 */
+TGraphErrors* GetFlow01_Rap02(Int_t n)
+{
+  // private communication 19.04.11, Raimond / Ante
+
+  if (n == 2)
+  {
+    //  v2{SP}(pt) for 0-1%, rapidity gap = 0.2:
+    const Int_t nPointsSP_0001ALICE_v2_etaGap02 = 18;
+    Double_t xSP_0001ALICE_v2_etaGap02[nPointsSP_0001ALICE_v2_etaGap02] = {0.250000,0.350000,0.450000,0.550000,0.650000,
+    0.750000,0.850000,0.950000,1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+    Double_t ySP_0001ALICE_v2_etaGap02[nPointsSP_0001ALICE_v2_etaGap02] = {0.009235,0.014105,0.017274,0.018245,0.023190,
+    0.024871,0.028216,0.031506,0.034179,0.035337,0.039836,0.045261,0.043393,0.050693,0.056469,0.055247,0.049718,0.052748};
+    Double_t xErrSP_0001ALICE_v2_etaGap02[nPointsSP_0001ALICE_v2_etaGap02] = {0.};
+    Double_t yErrSP_0001ALICE_v2_etaGap02[nPointsSP_0001ALICE_v2_etaGap02] = {0.000515,0.000504,0.000532,0.000585,0.000641,
+    0.000709,0.000788,0.000876,0.000723,0.000891,0.001098,0.001354,0.001671,0.001485,0.002463,0.004038,0.006441,0.008091};
+    TGraphErrors *SP_0001ALICE_v2_etaGap02 = new TGraphErrors(nPointsSP_0001ALICE_v2_etaGap02,xSP_0001ALICE_v2_etaGap02,
+                                                             ySP_0001ALICE_v2_etaGap02,xErrSP_0001ALICE_v2_etaGap02,yErrSP_0001ALICE_v2_etaGap02);
+    SP_0001ALICE_v2_etaGap02->SetMarkerStyle(kFullCircle);
+    SP_0001ALICE_v2_etaGap02->SetMarkerColor(kBlue);  
+    SP_0001ALICE_v2_etaGap02->SetFillStyle(1001);
+    SP_0001ALICE_v2_etaGap02->SetFillColor(kBlue-10);  
+    
+    return SP_0001ALICE_v2_etaGap02;
+  }
+  
+  if (n == 3)
+  {
+    //  v3{SP}(pt) for 0-1%, rapidity gap = 0.2:
+    const Int_t nPointsSP_0001ALICE_v3_etaGap02 = 18;
+    Double_t xSP_0001ALICE_v3_etaGap02[nPointsSP_0001ALICE_v3_etaGap02] = {0.250000,0.350000,0.450000,0.550000,0.650000,
+    0.750000,0.850000,0.950000,1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+    Double_t ySP_0001ALICE_v3_etaGap02[nPointsSP_0001ALICE_v3_etaGap02] = {0.005688,0.007222,0.010305,0.013795,0.016077,
+    0.018693,0.022310,0.026991,0.030162,0.035119,0.043097,0.048201,0.058249,0.063273,0.079233,0.083465,0.087807,0.069577};
+    Double_t xErrSP_0001ALICE_v3_etaGap02[nPointsSP_0001ALICE_v3_etaGap02] = {0.};
+    Double_t yErrSP_0001ALICE_v3_etaGap02[nPointsSP_0001ALICE_v3_etaGap02] = {0.000585,0.000582,0.000614,0.000667,0.000734,
+    0.000811,0.000898,0.000989,0.000817,0.001000,0.001234,0.001517,0.001874,0.001669,0.002765,0.004528,0.007202,0.009066};
+    TGraphErrors *SP_0001ALICE_v3_etaGap02 = new TGraphErrors(nPointsSP_0001ALICE_v3_etaGap02,xSP_0001ALICE_v3_etaGap02,ySP_0001ALICE_v3_etaGap02,
+                                                             xErrSP_0001ALICE_v3_etaGap02,yErrSP_0001ALICE_v3_etaGap02);
+    SP_0001ALICE_v3_etaGap02->SetMarkerStyle(kFullTriangleUp);
+    SP_0001ALICE_v3_etaGap02->SetMarkerSize(1.4);  
+    SP_0001ALICE_v3_etaGap02->SetMarkerColor(kGreen+2);
+    SP_0001ALICE_v3_etaGap02->SetFillStyle(1001);
+    SP_0001ALICE_v3_etaGap02->SetFillColor(kGreen-10);     
+    
+    return SP_0001ALICE_v3_etaGap02;
+  }
+   
+  if (n == 4)
+  {
+    //  v4{SP}(pt) for 0-1%, rapidity gap = 0.2:
+    const Int_t nPointsSP_0001ALICE_v4_etaGap02 = 18;
+    Double_t xSP_0001ALICE_v4_etaGap02[nPointsSP_0001ALICE_v4_etaGap02] = {0.250000,0.350000,0.450000,0.550000,0.650000,
+    0.750000,0.850000,0.950000,1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+    Double_t ySP_0001ALICE_v4_etaGap02[nPointsSP_0001ALICE_v4_etaGap02] = {0.001189,0.003540,0.004682,0.004210,0.007032,
+    0.008627,0.010226,0.013671,0.016214,0.020054,0.023878,0.033939,0.033693,0.040006,0.055612,0.066287,0.074857,0.078751};
+    Double_t xErrSP_0001ALICE_v4_etaGap02[nPointsSP_0001ALICE_v4_etaGap02] = {0.};
+    Double_t yErrSP_0001ALICE_v4_etaGap02[nPointsSP_0001ALICE_v4_etaGap02] = {0.001035,0.001017,0.001081,0.001187,0.001299,
+    0.001432,0.001590,0.001757,0.001443,0.001769,0.002175,0.002674,0.003296,0.002924,0.004844,0.007921,0.012524,0.015771};
+    TGraphErrors *SP_0001ALICE_v4_etaGap02 = new TGraphErrors(nPointsSP_0001ALICE_v4_etaGap02,xSP_0001ALICE_v4_etaGap02,ySP_0001ALICE_v4_etaGap02,
+                                                             xErrSP_0001ALICE_v4_etaGap02,yErrSP_0001ALICE_v4_etaGap02);
+    SP_0001ALICE_v4_etaGap02->SetMarkerStyle(kFullSquare);
+    SP_0001ALICE_v4_etaGap02->SetMarkerColor(kRed);
+    SP_0001ALICE_v4_etaGap02->SetFillStyle(1001);
+    SP_0001ALICE_v4_etaGap02->SetFillColor(kRed-10);  
+    
+    return SP_0001ALICE_v4_etaGap02;
+  }
+  
+  return 0;
+}
+
+TGraphErrors* GetFlow01_Rap10(Int_t n)
+{
+  if (n == 2)
+  {
+    //  v2{SP}(pt) for 0-1%, rapidity gap = 1.0:
+    const Int_t nPointsSP_0001ALICE_v2_etaGap10 = 21;
+    Double_t xSP_0001ALICE_v2_etaGap10[nPointsSP_0001ALICE_v2_etaGap10] = {0.250000,0.350000,0.450000,0.550000,0.650000,
+    0.750000,0.850000,0.950000,1.100000,1.300000,1.500000,1.700000,1.900000,2.100000,2.300000,2.500000,2.700000,2.900000,
+    3.250000,3.750000,4.500000};
+    Double_t ySP_0001ALICE_v2_etaGap10[nPointsSP_0001ALICE_v2_etaGap10] = {0.009129,0.013461,0.017567,0.018041,0.020384,
+    0.023780,0.021647,0.029543,0.028912,0.029464,0.037016,0.044131,0.043135,0.047286,0.051983,0.049311,0.050472,0.046569,
+    0.036905,0.054836,0.030527};
+    Double_t xErrSP_0001ALICE_v2_etaGap10[nPointsSP_0001ALICE_v2_etaGap10] = {0.};
+    Double_t yErrSP_0001ALICE_v2_etaGap10[nPointsSP_0001ALICE_v2_etaGap10] = {0.001179,0.001152,0.001219,0.001339,0.001480,
+    0.001644,0.001831,0.002016,0.001662,0.002033,0.002497,0.003056,0.003777,0.004645,0.005713,0.007069,0.008540,0.010447,
+    0.009145,0.014749,0.018698};
+    TGraphErrors *SP_0001ALICE_v2_etaGap10 = new TGraphErrors(nPointsSP_0001ALICE_v2_etaGap10,xSP_0001ALICE_v2_etaGap10,
+                                                             ySP_0001ALICE_v2_etaGap10,xErrSP_0001ALICE_v2_etaGap10,yErrSP_0001ALICE_v2_etaGap10);
+    SP_0001ALICE_v2_etaGap10->SetMarkerStyle(kOpenCircle);
+    SP_0001ALICE_v2_etaGap10->SetMarkerColor(kBlue);  
+    SP_0001ALICE_v2_etaGap10->SetFillStyle(1001);
+    SP_0001ALICE_v2_etaGap10->SetFillColor(kBlue-10);  
+    
+    return SP_0001ALICE_v2_etaGap10;
+  }
+
+  if (n == 3)
+  {
+    //  v3{SP}(pt) for 0-1%, rapidity gap = 1.0:
+    const Int_t nPointsSP_0001ALICE_v3_etaGap10 = 18;
+    Double_t xSP_0001ALICE_v3_etaGap10[nPointsSP_0001ALICE_v3_etaGap10] = {0.250000,0.350000,0.450000,0.550000,0.650000,
+    0.750000,0.850000,0.950000,1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000};
+    Double_t ySP_0001ALICE_v3_etaGap10[nPointsSP_0001ALICE_v3_etaGap10] = {0.006373,0.008403,0.010848,0.011505,0.016728,
+    0.018519,0.020163,0.027119,0.029315,0.036832,0.040974,0.043287,0.054395,0.060676,0.081763,0.074333,0.096016,0.074909};
+    Double_t xErrSP_0001ALICE_v3_etaGap10[nPointsSP_0001ALICE_v3_etaGap10] = {0.};
+    Double_t yErrSP_0001ALICE_v3_etaGap10[nPointsSP_0001ALICE_v3_etaGap10] = {0.001286,0.001269,0.001346,0.001474,0.001620,
+    0.001796,0.001991,0.002187,0.001803,0.002203,0.002697,0.003316,0.004078,0.003640,0.006050,0.009873,0.015824,0.020174};
+    TGraphErrors *SP_0001ALICE_v3_etaGap10 = new TGraphErrors(nPointsSP_0001ALICE_v3_etaGap10,xSP_0001ALICE_v3_etaGap10,ySP_0001ALICE_v3_etaGap10,
+                                                             xErrSP_0001ALICE_v3_etaGap10,yErrSP_0001ALICE_v3_etaGap10);
+    SP_0001ALICE_v3_etaGap10->SetMarkerStyle(kOpenTriangleUp);
+    SP_0001ALICE_v3_etaGap10->SetMarkerSize(1.2);  
+    SP_0001ALICE_v3_etaGap10->SetMarkerColor(kGreen+2);
+    SP_0001ALICE_v3_etaGap10->SetFillStyle(1001);
+    SP_0001ALICE_v3_etaGap10->SetFillColor(kGreen-10);     
+    
+    return SP_0001ALICE_v3_etaGap10;
+  }
+
+  if (n == 4)
+  {
+    //  v4{SP}(pt) for 0-1%, rapidity gap = 1.0:
+    const Int_t nPointsSP_0001ALICE_v4_etaGap10 = 11;
+    Double_t xSP_0001ALICE_v4_etaGap10[nPointsSP_0001ALICE_v4_etaGap10] = {0.300000,0.500000,0.700000,0.900000,1.200000,1.600000,2.000000,2.400000,2.800000,3.500000,4.500000};
+    Double_t ySP_0001ALICE_v4_etaGap10[nPointsSP_0001ALICE_v4_etaGap10] = {-0.000458,0.006444,0.005490,0.010870,0.018866,0.024370,0.029703,0.052505,0.060334,0.048189,0.128184};
+    Double_t xErrSP_0001ALICE_v4_etaGap10[nPointsSP_0001ALICE_v4_etaGap10] = {0.};
+    Double_t yErrSP_0001ALICE_v4_etaGap10[nPointsSP_0001ALICE_v4_etaGap10] = {0.001901,0.002012,0.002477,0.003014,0.002852,0.004297,0.006491,0.009846,0.014623,0.017120,0.040568};
+    TGraphErrors *SP_0001ALICE_v4_etaGap10 = new TGraphErrors(nPointsSP_0001ALICE_v4_etaGap10,xSP_0001ALICE_v4_etaGap10,ySP_0001ALICE_v4_etaGap10,
+                                                             xErrSP_0001ALICE_v4_etaGap10,yErrSP_0001ALICE_v4_etaGap10);
+    SP_0001ALICE_v4_etaGap10->SetMarkerStyle(kOpenSquare);
+    SP_0001ALICE_v4_etaGap10->SetMarkerColor(kRed);
+    SP_0001ALICE_v4_etaGap10->SetFillStyle(1001);
+    SP_0001ALICE_v4_etaGap10->SetFillColor(kRed-10);  
+    
+    return SP_0001ALICE_v4_etaGap10;
+  }
+}
+
+TGraphErrors* GetFlow02_Rap10(Int_t n)
+{
+  // private communication 20.04.11, Ante B. / Raimond
+
+  if (n == 2)
+  {
+     //  v2{SP}(pt) for 00-02%, eta gap = 1.0:
+    const Int_t nPointsSP_0002_v2_etaGap10 = 15;
+    Double_t xSP_0002_v2_etaGap10[nPointsSP_0002_v2_etaGap10] = {0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,
+    0.900000,1.100000,1.350000,1.650000,1.950000,2.250000,2.700000,3.500000,4.500000};
+    Double_t ySP_0002_v2_etaGap10[nPointsSP_0002_v2_etaGap10] = {0.010171,0.013190,0.017342,0.020629,0.022617,0.026549,
+    0.027423,0.032261,0.037467,0.041001,0.045763,0.049327,0.049688,0.051480,0.038527};
+    Double_t xErrSP_0002_v2_etaGap10[nPointsSP_0002_v2_etaGap10] = {0.};
+    Double_t yErrSP_0002_v2_etaGap10[nPointsSP_0002_v2_etaGap10] = {0.000600,0.000590,0.000625,0.000683,0.000757,0.000839,
+    0.000692,0.000848,0.000888,0.001209,0.001653,0.002252,0.002465,0.003968,0.009391};
+    TGraphErrors *SP_0002_v2_etaGap10 = new TGraphErrors(nPointsSP_0002_v2_etaGap10,xSP_0002_v2_etaGap10,ySP_0002_v2_etaGap10,
+                                                  xErrSP_0002_v2_etaGap10,yErrSP_0002_v2_etaGap10);
+                                                 
+    return SP_0002_v2_etaGap10;
+  }
+  
+  if (n == 3)
+  {
+    const Int_t nPointsSP_0002_v3_etaGap10 = 15;
+    Double_t xSP_0002_v3_etaGap10[nPointsSP_0002_v3_etaGap10] = {0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,
+    0.900000,1.100000,1.350000,1.650000,1.950000,2.250000,2.700000,3.500000,4.500000};
+    Double_t ySP_0002_v3_etaGap10[nPointsSP_0002_v3_etaGap10] = {0.006592,0.007286,0.012180,0.012242,0.017416,0.018393,
+    0.024716,0.030980,0.037703,0.046558,0.051285,0.064613,0.074831,0.077093,0.082442};
+    Double_t xErrSP_0002_v3_etaGap10[nPointsSP_0002_v3_etaGap10] = {0.};
+    Double_t yErrSP_0002_v3_etaGap10[nPointsSP_0002_v3_etaGap10] = {0.000682,0.000676,0.000713,0.000782,0.000860,0.000953,
+    0.000782,0.000957,0.001002,0.001361,0.001862,0.002541,0.002767,0.004466,0.010586};
+    TGraphErrors *SP_0002_v3_etaGap10 = new TGraphErrors(nPointsSP_0002_v3_etaGap10,xSP_0002_v3_etaGap10,ySP_0002_v3_etaGap10,
+                                                         xErrSP_0002_v3_etaGap10,yErrSP_0002_v3_etaGap10);    
+                                                         
+    return SP_0002_v3_etaGap10;
+  }
+  
+  if (n == 4)
+  {
+    const Int_t nPointsSP_0002_v4_etaGap10 = 15;
+    Double_t xSP_0002_v4_etaGap10[nPointsSP_0002_v4_etaGap10] = {0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,
+    0.900000,1.100000,1.350000,1.650000,1.950000,2.250000,2.700000,3.500000,4.500000};
+    Double_t ySP_0002_v4_etaGap10[nPointsSP_0002_v4_etaGap10] = {-0.000533,0.001167,0.002081,0.005218,0.006826,0.008440,
+    0.013009,0.014812,0.017125,0.030106,0.038279,0.050488,0.067640,0.071637,0.084239};
+    Double_t xErrSP_0002_v4_etaGap10[nPointsSP_0002_v4_etaGap10] = {0.};
+    Double_t yErrSP_0002_v4_etaGap10[nPointsSP_0002_v4_etaGap10] = {0.001427,0.001398,0.001482,0.001594,0.001758,0.001945,
+    0.001593,0.001951,0.002046,0.002787,0.003802,0.005182,0.005663,0.009064,0.021449};
+    TGraphErrors *SP_0002_v4_etaGap10 = new TGraphErrors(nPointsSP_0002_v4_etaGap10,xSP_0002_v4_etaGap10,ySP_0002_v4_etaGap10,
+                                                      xErrSP_0002_v4_etaGap10,yErrSP_0002_v4_etaGap10);
+    return SP_0002_v4_etaGap10;
+  }
+}
 
 TGraphErrors* GetFlow02(Int_t n)
 {
@@ -1948,8 +2158,10 @@ TGraphErrors* GetFlow4050(Int_t n)
 
 Float_t CalculateFlow(TH1* ptDist, Float_t ptMin, Float_t ptMax, Int_t n, Int_t centralityBegin, Int_t centralityEnd)
 {
-  if (centralityBegin == 0 && centralityEnd == 2)
-    flow = GetFlow02(n);
+  if (centralityBegin == 0 && centralityEnd == 1)
+    flow = GetFlow01_Rap10(n);
+  else if (centralityBegin == 0 && centralityEnd == 2)
+    flow = GetFlow02_Rap10(n);
   else if (centralityBegin == 0 && centralityEnd == 5)
     flow = GetFlow05(n);
   else if (centralityBegin == 5 && centralityEnd == 10)
@@ -2155,7 +2367,8 @@ void GraphShiftX(TGraphErrors* graph, Float_t offset)
     graph->GetX()[i] += offset;
 }
 
-const Float_t kPythiaScalingFactor = 0.935;
+// Float_t kPythiaScalingFactor = 0.935;
+Float_t kPythiaScalingFactor = 1;
 
 void DrawYields(const char* fileName = "yields.root")
 {
@@ -3037,8 +3250,10 @@ void CompareIAAICP(const char* fileName1, const char* fileName2, Int_t nominator
   //Int_t caseId = 12;
 
   nearSide1 = GetRatio(fileName1, nominatorBin, denominatorBin, j, caseId, 0);
-  nearSide2 = GetRatio(fileName2, nominatorBin, denominatorBin, j, caseId, 0);
   awaySide1 = GetRatio(fileName1, nominatorBin, denominatorBin, j, caseId, 1);
+  
+  kPythiaScalingFactor = 1;
+  nearSide2 = GetRatio(fileName2, nominatorBin, denominatorBin, j, caseId, 0);
   awaySide2 = GetRatio(fileName2, nominatorBin, denominatorBin, j, caseId, 1);
   
 /*  ScaleGraph(nearSide2, 1.33);
@@ -3246,6 +3461,9 @@ void IAA(const char* fileName, Int_t iaa)
   // 0 = IAA LHC1
   // 1 = IAA RHIC
   // 2 = ICP
+  
+  if (kPythiaScalingFactor != 1)
+    Printf("Using reference data scaling factor: %f", kPythiaScalingFactor);
 
   style();
   
@@ -3463,7 +3681,8 @@ void IAA(const char* fileName, Int_t iaa)
       latex->Draw();
     }
   
-    latex = new TLatex(0.5+xC, 0.90, "ALICE preliminary");
+//     latex = new TLatex(0.5+xC, 0.90, "ALICE preliminary");
+    latex = new TLatex(0.5+xC, 0.90, "-- work in progress --");
     latex->SetTextSize(0.04);
     latex->SetNDC();
     latex->Draw();
@@ -3594,14 +3813,16 @@ Float_t ExtractYields(TH1* hist, Int_t trigId, Int_t centralityId, Float_t ptA,
   Int_t binBegin = hist->FindBin(-width);
   Int_t binEnd = hist->FindBin(width);
   Double_t nearSideError = 0;
-  Float_t nearSideYield = hist->IntegralAndError(binBegin, binEnd, nearSideError) / (binEnd - binBegin + 1);
-  nearSideError /= (binEnd - binBegin + 1);
-  
+  Float_t nearSideYield = hist->IntegralAndError(binBegin, binEnd, nearSideError);
+  nearSideYield *= hist->GetXaxis()->GetBinWidth(1);
+  nearSideError *= hist->GetXaxis()->GetBinWidth(1);
+
   Int_t binBegin = hist->FindBin(TMath::Pi() - width);
   Int_t binEnd = hist->FindBin(TMath::Pi() + width);
   Double_t awaySideError = 0;
-  Float_t awaySideYield = hist->IntegralAndError(binBegin, binEnd, awaySideError) / (binEnd - binBegin + 1);
-  awaySideError /= (binEnd - binBegin + 1);
+  Float_t awaySideYield = hist->IntegralAndError(binBegin, binEnd, awaySideError);
+  awaySideYield *= hist->GetXaxis()->GetBinWidth(1);
+  awaySideError *= hist->GetXaxis()->GetBinWidth(1);
   
   Printf("%d %d %f %d: %f+-%f %f+-%f (%f)", trigId, centralityId, ptA, caseId, nearSideYield, nearSideError, awaySideYield, awaySideError, normUnc);
   
@@ -4067,13 +4288,120 @@ void RemoveBaseLine(TH1* hist)
     hist->SetBinContent(i, hist->GetBinContent(i) - zyam);
 }
  
+TH2* cacheMixed = 0;
 Int_t gHistCount = 0;
-void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t& v2, Int_t step, Int_t centralityBegin, Int_t centralityEnd, Float_t ptBegin, Float_t ptEnd, Int_t twoD = 0)
+void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t* v2, Int_t step, Int_t centralityBegin, Int_t centralityEnd, Float_t ptBegin, Float_t ptEnd, Int_t twoD = 0, Bool_t equivMixedBin = kFALSE, Float_t* vn = 0, Bool_t scaleToPairs = kTRUE)
 {
   h = (AliUEHistograms*) hVoid;
   hMixed = (AliUEHistograms*) hMixedVoid;
 
-  *hist = h->GetUEHist(2)->GetUEHist(step, 0, ptBegin, ptEnd, h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(0.01 + centralityBegin), h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(-0.01 + centralityEnd), twoD);
+  Int_t centralityBeginBin = 0;
+  Int_t centralityEndBin = -1;
+  
+  if (centralityEnd >= centralityBegin)
+  {
+    centralityBeginBin = h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(0.01 + centralityBegin);
+    centralityEndBin = h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(-0.01 + centralityEnd);
+  }
+  
+  // 2d same and mixed event
+  TH2* sameTwoD  = h->GetUEHist(2)->GetUEHist(step, 0, ptBegin, ptEnd, centralityBeginBin, centralityEndBin, 1, kFALSE);
+  
+  if (hMixed)
+  {
+    if (!equivMixedBin)
+    {
+      // No centrality, nor pT associated dep of the mixed event observed. Use a larger range to get more statistics
+      
+      // TODO HACK. mixed event is not propagated to step0
+      Int_t stepMixed = 6;
+      if (!cacheMixed)
+      {
+       hMixed->SetPtRange(1.0, 10);
+       cacheMixed = (TH2*) hMixed->GetUEHist(2)->GetUEHist(stepMixed, 0, 1.0, 10.0, 1, 15, 1, kFALSE);
+      }
+      
+      TH2* mixedTwoD = cacheMixed;
+    }
+    else
+    {
+      // use same bin for mixing
+      
+      Int_t stepMixed = 6;
+      TH2* mixedTwoD = (TH2*) hMixed->GetUEHist(2)->GetUEHist(stepMixed, 0, ptBegin, ptEnd, centralityBeginBin, centralityEndBin, 1, kFALSE);
+    }
+    
+    if (0)
+    {
+      // asssume flat in dphi, gain in statistics
+      
+      TH1* histMixedproj = mixedTwoD->ProjectionY();
+      histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
+      
+      for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
+       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
+         mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
+    }
+    
+    // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
+    if (scaleToPairs)
+    {
+      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
+      mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
+    }
+    else
+      Double_t mixedNorm = mixedTwoD->Integral() / sameTwoD->Integral();
+    
+    // divide and scale
+    sameTwoD->Divide(mixedTwoD);
+    sameTwoD->Scale(mixedNorm);
+    
+/*    new TCanvas;
+    sameTwoD->Draw("SURF1");
+    dfdsafd;*/
+  }
+  
+  // extract dphi distribution if requested
+  if (twoD == 1)
+  {
+    *hist = sameTwoD;
+  }
+  else
+  {
+    Float_t etaLimit = 0.8;
+
+    if (twoD == 0 || twoD == 10)
+    {
+      Int_t etaBegin = 1;
+      Int_t etaEnd = sameTwoD->GetNbinsY();
+      
+      if (twoD == 10)
+      {
+       etaBegin = sameTwoD->GetYaxis()->FindBin(-etaLimit + 0.01);
+       etaEnd   = sameTwoD->GetYaxis()->FindBin(etaLimit - 0.01);
+      }
+
+      *hist = sameTwoD->ProjectionX("proj", etaBegin, etaEnd);
+      
+      if (!scaleToPairs)
+       (*hist)->Scale(1.0 / (etaEnd - etaBegin + 1));
+    }
+    else if (twoD == 11)
+    {
+      *hist = sameTwoD->ProjectionX("proj", sameTwoD->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), sameTwoD->GetYaxis()->FindBin(-etaLimit - 0.01));
+      Int_t etaBins = sameTwoD->GetYaxis()->FindBin(-etaLimit - 0.01) - sameTwoD->GetYaxis()->FindBin(-etaLimit * 2 + 0.01) + 1;
+
+      TH1D* tracksTmp = sameTwoD->ProjectionX("proj2", sameTwoD->GetYaxis()->FindBin(etaLimit + 0.01), sameTwoD->GetYaxis()->FindBin(2 * etaLimit - 0.01));
+      etaBins += sameTwoD->GetYaxis()->FindBin(2 * etaLimit - 0.01) - sameTwoD->GetYaxis()->FindBin(etaLimit + 0.01) + 1;
+
+      (*hist)->Add(tracksTmp);
+      
+      if (!scaleToPairs)
+       (*hist)->Scale(1.0 / etaBins);
+    }
+  }
+  
+  //*hist = h->GetUEHist(2)->GetUEHist(step, 0, ptBegin, ptEnd, h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(0.01 + centralityBegin), h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(-0.01 + centralityEnd), twoD);
   
   *hist = (TH1*) (*hist)->Clone(Form("GetDistAndFlow_%d", gHistCount++));
   
@@ -4087,7 +4415,7 @@ void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t& v2, Int_
   newTitle.Form("%s - %s - %d-%d%%", str.Data(), str2.Data(), centralityBegin, centralityEnd);
   (*hist)->SetTitle(newTitle);
   
-  if (hMixed)
+  if (0 && hMixed)
   {
     histMixed = hMixed->GetUEHist(2)->GetUEHist(step, 0, ptBegin, ptEnd, hMixed->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(0.01 + centralityBegin), hMixed->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(-0.01 + centralityEnd));
     
@@ -4101,21 +4429,32 @@ void GetDistAndFlow(void* hVoid, void* hMixedVoid, TH1** hist, Float_t& v2, Int_
     //(*hist)->DrawCopy("SAME")->SetLineColor(4);
   }
   
-  // calculate v2trigger
-  h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->SetRangeUser(0.01 + centralityBegin, -0.01 + centralityEnd);
-  ptDist = h->GetUEHist(2)->GetEventHist()->Project(step, 0);
-  Float_t v2Trig = CalculateFlow(ptDist, ptBegin, ptEnd, 2, centralityBegin, centralityEnd);
-  delete ptDist;
-  
-  // calculate v2 assoc
-  cont = h->GetUEHist(2)->GetTrackHist(0);
-  h->GetUEHist(2)->GetTrackHist(0)->GetGrid(step)->GetGrid()->GetAxis(3)->SetRangeUser(0.01 + centralityBegin, -0.01 + centralityEnd);
-  h->GetUEHist(2)->GetTrackHist(0)->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(ptBegin, ptEnd);
-  ptDist = h->GetUEHist(2)->GetTrackHist(0)->Project(step, 1);
-  Float_t v2Assoc = CalculateFlow(ptDist, gpTMin, gpTMax, 2, centralityBegin, centralityEnd);
-  delete ptDist;
-  
-  v2 = v2Trig * v2Assoc;
+  if (v2 || vn)
+  {
+    // calculate v2trigger
+    h->GetUEHist(2)->GetEventHist()->GetGrid(step)->GetGrid()->GetAxis(1)->SetRangeUser(0.01 + centralityBegin, -0.01 + centralityEnd);
+    ptDist = h->GetUEHist(2)->GetEventHist()->Project(step, 0);
+    Float_t vTrig[4];
+    for (Int_t i=2; i<=((vn) ? 4 : 2); i++)
+      vTrig[i-1] = CalculateFlow(ptDist, ptBegin, ptEnd, i, centralityBegin, centralityEnd);
+    delete ptDist;
+    
+    // calculate v2 assoc
+    cont = h->GetUEHist(2)->GetTrackHist(0);
+    h->GetUEHist(2)->GetTrackHist(0)->GetGrid(step)->GetGrid()->GetAxis(3)->SetRangeUser(0.01 + centralityBegin, -0.01 + centralityEnd);
+    h->GetUEHist(2)->GetTrackHist(0)->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(ptBegin, ptEnd);
+    ptDist = h->GetUEHist(2)->GetTrackHist(0)->Project(step, 1);
+    Float_t vAssoc[4];
+    for (Int_t i=2; i<=((vn) ? 4 : 2); i++)
+      vAssoc[i-1] = CalculateFlow(ptDist, gpTMin, gpTMax, i, centralityBegin, centralityEnd);
+    delete ptDist;
+  
+    if (v2)
+      *v2 = vTrig[2-1] * vAssoc[2-1];
+    if (vn)
+      for (Int_t i=2; i<=4; i++)
+       vn[i-1] = vTrig[i-1] * vAssoc[i-1];
+  }
 }
  
 void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Float_t yMax = 0.1, Int_t twoD = 0, Int_t centrBegin = 1, Int_t centrEnd = 1)
@@ -4125,7 +4464,7 @@ void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Flo
   Bool_t veryCentral = 0;
   Bool_t flowComparison = 0;
   Bool_t rhicOverlay = 0;
-  Bool_t highStatBinning = 1;
+  Bool_t highStatBinning = 0;
 
   file = TFile::Open("dphi_corr.root", "RECREATE");
   file->Close();
@@ -4183,7 +4522,13 @@ void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Flo
   AliUEHistograms* hMixed = 0;
   if (twoD)
     hMixed = (AliUEHistograms*) GetUEHistogram(fileName1, 0, kTRUE);
-    
+
+  if (veryCentral)
+  {
+    Printf("WARNING: Reading mixed event from preliminaries/corrected_110317.root");
+    hMixed = (AliUEHistograms*) GetUEHistogram("preliminaries/corrected_110317.root", 0, kTRUE);
+  }
+  
   AliUEHistograms* h2 = 0;
   if (!twoD)
     h2 = (AliUEHistograms*) GetUEHistogram(fileName2);
@@ -4272,11 +4617,19 @@ void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Flo
         
         if (veryCentral)
         {
-          Int_t step = 6;
-          TH1* hist1 = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 1, 2);  hist1Str = "0-2%";
-          TH1* hist2 = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 2, 3);  hist2Str = "2-3%";
-          TH1* hist2b = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 10, 10); hist2bStr = "30-40%";
-          TH1* hist3 = 0; // h2->GetUEHist(2)->GetUEHist(0, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01);
+          Int_t step = 0;
+          TH1* hist1 = 0;
+          TH1* hist2 = 0;
+          TH1* hist2b = 0;
+         TH1* hist3 = 0;
+         
+         GetDistAndFlow(h, hMixed, &hist1,  v2, step, 0,  2,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist1Str = "0-2%";
+//       GetDistAndFlow(h, hMixed, &hist2,  v2, step, 1,  3,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2Str = "1-3%";
+         GetDistAndFlow(h, hMixed, &hist2b,  v2+2, step, 30,  40,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2bStr = "30-40%";
+          
+         //TH1* hist1 = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 1, 2);  hist1Str = "0-2%";
+          //TH1* hist2 = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 2, 3);  hist2Str = "2-3%";
+          //TH1* hist2b = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 10, 10); hist2bStr = "30-40%";
         }
         else if (flowComparison)
         {
@@ -4294,14 +4647,14 @@ void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Flo
         }
         else
         {
-          Int_t step = 6;
+          Int_t step = 0;
           TH1* hist1 = 0;
           TH1* hist2 = 0;
           TH1* hist2b = 0;
-          GetDistAndFlow(h, hMixed, &hist1,  v2[0], step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist1Str = "0-5%";
-          GetDistAndFlow(h, hMixed, &hist2,  v2[1], step, 0, 20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2Str = "0-20%";
+          GetDistAndFlow(h, hMixed, &hist1,  v2, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist1Str = "0-5%";
+          GetDistAndFlow(h, hMixed, &hist2,  v2+1, step, 0, 20, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2Str = "0-20%";
 //           GetDistAndFlow(h, hMixed, &hist2b, v2[2], step, 60, 80, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2bStr = "60-80%";
-          GetDistAndFlow(h, hMixed, &hist2b, v2[2], step, 60, 90, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2bStr = "60-90%";
+          GetDistAndFlow(h, hMixed, &hist2b, v2+2, step, 60, 90, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01); hist2bStr = "60-90%";
           
           Printf("%f %f %f", v2[0], v2[1], v2[2]);
           
@@ -4309,9 +4662,15 @@ void PlotDeltaPhiDistributions(const char* fileName1, const char* fileName2, Flo
 //           TH1* hist2 = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 9, 12);   hist2Str = "20-60%";
 //           TH1* hist2b = h->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, 13, 15); hist2bStr = "60-90%";
           
-          step = 6;
-          TH1* hist3 = h2->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01);
-        }
+          step = 0;
+//           TH1* hist3Old = h2->GetUEHist(2)->GetUEHist(step, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01);
+         TH1* hist3 = 0;
+          GetDistAndFlow(h2, hMixed, &hist3,  0, step, 0,  -1,  leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01);
+         
+/*       new TCanvas;
+         hist3->Draw();
+         hist3Old->DrawCopy("SAME")->SetLineColor(2); */
+       }
         
         /*
         RemoveBaseLine(hist1);
@@ -5377,8 +5736,8 @@ void DeltaPhiBaseLinePreliminary(const char* fileName)
   TH1* hist1 = 0;
   TH1* hist2b = 0;
   
-  GetDistAndFlow(h, 0, &hist1,  v2[0], step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
-  GetDistAndFlow(h, 0, &hist2b, v2[2], step, 60, 90, leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
+  GetDistAndFlow(h, 0, &hist1,  v2, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
+  GetDistAndFlow(h, 0, &hist2b, v2+2, step, 60, 90, leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
   
   hist1->SetLineWidth(2);
   hist2b->SetLineWidth(2);
@@ -5472,17 +5831,23 @@ void DeltaPhiRidgePreliminary(const char* fileName)
   
   loadlibs();
 
-  Float_t leadingPtArr[] = { 6.0, 8.0, 10.0, 15.0, 15.0 };
-  Float_t assocPtArr[] =     { 0.5, 1.5, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0 };
-  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+  Float_t leadingPtArr[] = { 2.0, 2.0, 3.0, 4.0, 10.0, 20.0, 40.0 };
+  Float_t assocPtArr[] =   { 1.0, 2.0, 3.0, 6.0, 10.0, 20.0, 40.0 };
+//   Float_t leadingPtArr[] = { 6.0, 8.0, 10.0, 15.0, 15.0 };
+//   Float_t assocPtArr[] =     { 0.5, 1.5, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0 };
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName, 0, kFALSE);
+  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
 
-  Int_t i = 1;
-  Int_t j = 3;
-  Int_t step = 0;
+//   Int_t i = 1;
+//   Int_t j = 3;
+  Int_t i = 0;
+  Int_t j = 0;
+  Int_t step = 6;
 
   gpTMin = assocPtArr[j] + 0.01;
   gpTMax = assocPtArr[j+1] - 0.01;
   SetupRanges(h);
+  SetupRanges(hMixed);
       
   TString str;
   str.Form("%.1f < p_{T,trig} < %.1f", leadingPtArr[i], leadingPtArr[i+2]);
@@ -5498,11 +5863,20 @@ void DeltaPhiRidgePreliminary(const char* fileName)
   TH1* hist1Ridge = 0;
   TH1* hist2bRidge = 0;
   
-  GetDistAndFlow(h, 0, &hist1,  v2[0], step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
+  TH1* hist2d = 0;
+  TH1* hist2dMixed = 0;
+  
+//   GetDistAndFlow(h, 0, &hist2d,  v2, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 1); 
+//   GetDistAndFlow(hMixed, 0, &hist2dMixed,  v2, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 1); 
+//   new TCanvas;
+//   //hist2d->Divide(hist2dMixed);
+//   hist2dMixed->Draw("SURF1");
+  
+  GetDistAndFlow(h, hMixed, &hist1,  0, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
 //   GetDistAndFlow(h, 0, &hist2b, v2[2], step, 60, 90, leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01); 
-  GetDistAndFlow(h, 0, &hist1Peak,  v2[0], step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 10); 
+  GetDistAndFlow(h, hMixed, &hist1Peak,  0, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 10); 
 
-  GetDistAndFlow(h, 0, &hist1Ridge,  v2[0], step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 11); 
+  GetDistAndFlow(h, hMixed, &hist1Ridge,  0, step, 0,  5,  leadingPtArr[i] + 0.01, leadingPtArr[i+2] - 0.01, 11); 
 
   // TODO normalize
 
@@ -5517,9 +5891,6 @@ void DeltaPhiRidgePreliminary(const char* fileName)
   hist1->SetTitle("");
   hist1->DrawCopy();
   
-/*  hist1Peak->Scale(1./0.75);
-  hist1Ridge->Scale(1./0.25);*/
-  
   hist1Peak->SetLineColor(2);
   hist1Peak->DrawCopy("SAME");
   
@@ -5699,7 +6070,7 @@ void CompareDeltaPhi(const char* fileName1, const char* fileName2, Int_t central
     }
 }
   
-TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption, Bool_t plotContributions, Int_t twoD)
+TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption, Bool_t plotContributions, Int_t twoD, TF1** flowFuncP = 0)
 {
   if (0)
   {
@@ -5729,18 +6100,22 @@ TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption,
     loadlibs();
     
     AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+    AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
     
     Float_t leadingPtArr[] = { 2.0, 3.0, 4.0, 10.0, 20.0, 40.0 };
     Float_t assocPtArr[] =   { 1.0, 2.0, 3.0, 6.0, 10.0, 20.0, 40.0 };
+//     Float_t assocPtArr[] =   { 2.0, 3.0, 6.0, 10.0, 20.0, 40.0 };
   
     Int_t i = centrality % 2;
     Int_t step = 6;
     
+    Int_t j = 0;
+    
     gpTMin = assocPtArr[i] + 0.01;
     gpTMax = assocPtArr[i+1] - 0.01;
     
     Int_t centralityBegin = 0;
-    Int_t centralityEnd = 2;
+    Int_t centralityEnd = 1;
     
     if (centrality >= 2)
     {
@@ -5749,13 +6124,54 @@ TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption,
     }
       
     SetupRanges(h);
+    SetupRanges(hMixed);
     
-    Float_t v2[1];
     TH1* hist = 0;
-    GetDistAndFlow(h, 0, &hist,  v2[0], step, centralityBegin, centralityEnd, leadingPtArr[i] + 0.01, leadingPtArr[i+1] - 0.01, twoD); 
+    Float_t vn[4];
+    vn[0] = 0;
+    vn[1] = 0;
+    vn[2] = 0;
+    vn[3] = 0;
+    
+    Bool_t scaleToPairs = 1;
+    
+    GetDistAndFlow(h, hMixed, &hist, 0, step, centralityBegin, centralityEnd, leadingPtArr[j] + 0.01, leadingPtArr[j+1] - 0.01, twoD, kTRUE, vn, scaleToPairs); 
+    if (scaleToPairs && twoD == 11)
+      hist->Scale(2);  
+    
+    Printf("%f %f %f", vn[2], vn[3], vn[4]);
+    
+    // mirror delta phi to improve stats
+    if (1)
+    {
+      for (Int_t bin=1; bin<=hist->GetNbinsX(); bin++)
+      {
+       if (hist->GetBinCenter(bin) < 0 || hist->GetBinCenter(bin) > TMath::Pi())
+       {
+         if (hist->GetBinCenter(bin) < 0)
+           Int_t bin2 = hist->FindBin(-1.0 * hist->GetBinCenter(bin));
+         else
+           Int_t bin2 = hist->FindBin(TMath::TwoPi() - hist->GetBinCenter(bin));
+         
+         Float_t combValue = hist->GetBinContent(bin) / hist->GetBinError(bin) / hist->GetBinError(bin) + hist->GetBinContent(bin2) / hist->GetBinError(bin2) / hist->GetBinError(bin2);
+         
+         Float_t weight = 1. / hist->GetBinError(bin) / hist->GetBinError(bin) + 1. / hist->GetBinError(bin2) / hist->GetBinError(bin2);
+         combValue /= weight;
+         
+         Float_t combError = TMath::Sqrt(1.0 / weight);
+         
+         hist->SetBinContent(bin2, combValue);
+         hist->SetBinError(bin2, combError);
+
+         hist->SetBinContent(bin, combValue);
+         hist->SetBinError(bin, combError);
+         
+         Printf("%d %d %f %f", bin, bin2, combValue, combError);
+       }
+      }
+    }
   }
   
-  
   TString str(hist->GetTitle());
   str.ReplaceAll(" - ", "#");
   tokens = str.Tokenize("#");
@@ -5764,90 +6180,24 @@ TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption,
   hist->SetLineColor(1);
   hist->SetYTitle("1/N_{trig} dN_{assoc}/d#Delta#phi");
   hist->SetXTitle("#Delta#phi (rad.)");
+  hist = (TH1*) hist->Clone();
+  if (!drawingOption)
+    return hist;
+  
   hist = hist->DrawCopy(drawingOption);
   
-  for (Int_t i=0; i<tokens->GetEntries(); i++)
+  for (Int_t i=0; i<tokens->GetEntries()-1; i++)
   {
-    latex = new TLatex(0.6, 0.8-i*0.05, tokens->At(i)->GetName());
+    if (centrality == 1)
+      latex = new TLatex(0.6, 0.9 - i*0.05 - gPad->GetTopMargin(), tokens->At(i)->GetName());
+    else
+      latex = new TLatex(0.6, 0.12 - i*0.05, tokens->At(i)->GetName());
+    
     latex->SetNDC();
     latex->SetTextSize(0.04);
     latex->Draw();
   }
   
-  // From CalculateFlow(...) above
-  // v_{2} = 0.042293 for 1.000000 < pT < 2.000000
-  // v_{2} = 0.060797 for 2.000000 < pT < 4.000000
-  // v_{3} = 0.040318 for 1.000000 < pT < 2.000000
-  // v_{3} = 0.077353 for 2.000000 < pT < 4.000000
-  // v_{4} = 0.020557 for 1.000000 < pT < 2.000000
-  // v_{4} = 0.056187 for 2.000000 < pT < 4.000000
-  //Double_t vn[] = { 0, 0.042293 * 0.060797, 0.040318 * 0.077353, 0.020557 * 0.056187 };
-  
-  switch (centrality)
-  {
-    case 0: // 0-2%, 1-2 times 2-3
-/*      v_{2} = 0.042240 for 1.000000 < pT < 2.000000
-      v_{2} = 0.060385 for 2.000000 < pT < 3.000000
-      v_{3} = 0.040230 for 1.000000 < pT < 2.000000
-      v_{3} = 0.074648 for 2.000000 < pT < 3.000000
-      v_{4} = 0.020487 for 1.000000 < pT < 2.000000
-      v_{4} = 0.053341 for 2.000000 < pT < 3.000000*/
-      Double_t vn[] = { 0, 0.042240 * 0.060385, 0.040230 * 0.074648, 0.020487 * 0.053341 };
-      break;
-  
-    case 1: // 0-2%, 2-3 times 3-4
-      // v_{2} = 0.058958 for 2.000000 < pT < 3.000000
-      // v_{2} = 0.063784 for 3.000000 < pT < 4.000000
-      // v_{3} = 0.072185 for 2.000000 < pT < 3.000000
-      // v_{3} = 0.096924 for 3.000000 < pT < 4.000000
-      // v_{4} = 0.049694 for 2.000000 < pT < 3.000000
-      // v_{4} = 0.076772 for 3.000000 < pT < 4.000000
-      Double_t vn[] = { 0, 0.058958 * 0.063784, 0.072185 * 0.096924, 0.049694 * 0.076772 };
-      break;
-      
-      
-//     case 1: 
-//    /* v_{2} = 0.113591 for 2.000000 < pT < 3.000000
-//       v_{2} = 0.130848 for 3.000000 < pT < 4.000000
-//       v_{3} = 0.085376 for 2.000000 < pT < 3.000000
-//       v_{3} = 0.109396 for 3.000000 < pT < 4.000000
-//       v_{4} = 0.056515 for 2.000000 < pT < 3.000000
-//       v_{4} = 0.087579 for 3.000000 < pT < 4.000000*/
-//       Double_t vn[] = { 0, 0.113591 * 0.130848, 0.085376 * 0.109396, 0.056515 * 0.087579 };
-//       break;
-//       
-//     case 2:
-//       /*v_{2} = 0.205376 for 2.000000 < pT < 3.000000
-//       v_{2} = 0.224529 for 3.000000 < pT < 4.000000
-//       v_{3} = 0.104202 for 2.000000 < pT < 3.000000
-//       v_{3} = 0.124636 for 3.000000 < pT < 4.000000
-//       v_{4} = 0.069736 for 2.000000 < pT < 3.000000
-//       v_{4} = 0.098462 for 3.000000 < pT < 4.000000*/
-//       Double_t vn[] = { 0, 0.205376 * 0.224529, 0.104202 * 0.124636, 0.069736 * 0.098462 };
-//       break;
-//       
-   case 2: // 30-40%, 1-2 times 2-3
-      /*v_{2} = 0.170185 for 1.000000 < pT < 2.000000
-      v_{2} = 0.236121 for 2.000000 < pT < 3.000000
-      v_{3} = 0.068512 for 1.000000 < pT < 2.000000
-      v_{3} = 0.115416 for 2.000000 < pT < 3.000000
-      v_{4} = 0.037889 for 1.000000 < pT < 2.000000
-      v_{4} = 0.079676 for 2.000000 < pT < 3.000000*/
-      Double_t vn[] = { 0, 0.170185 * 0.236121, 0.068512 * 0.115416, 0.037889 * 0.079676 };
-      break;
-   
-   case 3: // 30-40%, 2-3 times 3-4
-      /*v_{2} = 0.231893 for 2.000000 < pT < 3.000000
-      v_{2} = 0.243978 for 3.000000 < pT < 4.000000
-      v_{3} = 0.111616 for 2.000000 < pT < 3.000000
-      v_{3} = 0.131851 for 3.000000 < pT < 4.000000
-      v_{4} = 0.075781 for 2.000000 < pT < 3.000000
-      v_{4} = 0.107771 for 3.000000 < pT < 4.000000*/
-
-      Double_t vn[] = { 0, 0.231893 * 0.243978, 0.111616 * 0.131851, 0.075781 * 0.107771 };
-      break;
-  }
-  
   TF1* flowFunc = new TF1("flowFunc", "[0] * (1+2*[1]*cos(2*x)+2*[2]*cos(3*x)+2*[3]*cos(4*x))", -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   //flowFunc->SetLineWidth(1);
   flowFunc->SetLineColor(2);
@@ -5859,11 +6209,15 @@ TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption,
   //flowFunc->SetParameter(0, 85.4); flowFunc->DrawClone("SAME");
   //flowFunc->SetParameter(0, 85.6); flowFunc->DrawClone("SAME");
   
-  hist->Fit(flowFunc, "", "SAME", 1.8, 4.5);
+//   hist->Fit(flowFunc, "", "SAME", 1.8, 4.5);
+  hist->Fit(flowFunc, "N", "", -0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   hist->GetYaxis()->SetRangeUser(hist->GetMinimum() * 0.95, hist->GetMaximum() * 1.05);
   flowFunc->SetRange(-0.5 * TMath::Pi(), 1.5 * TMath::Pi());
   flowFunc->Draw("SAME");
   
+  if (flowFuncP)
+    *flowFuncP = flowFunc;
+  
   //hist->GetYaxis()->SetRangeUser(84, 89);
   //flowFunc->SetParameter(4, 1);
   
@@ -5895,33 +6249,118 @@ TH1* MACHCone(const char* fileName, Int_t centrality, const char* drawingOption,
   return hist;
 }
 
-void MACHConeAll(const char* fileName, Int_t twoD)
+void MACHConeAll(const char* fileName)
 {
   loadlibs();
 
   SetFlowStyle();
 
-  c = new TCanvas("c", "c", 800, 800);
-  c->Divide(2, 2);
+  c = new TCanvas("c", "c", 300, 800);
+  c->SetTopMargin(0);
+  c->SetLeftMargin(0);
+  c->SetRightMargin(0);
+  c->SetBottomMargin(0);
+  c->Divide(1, 2, 0, 0);
   
-  Float_t rangesMin[] = { 42, 5.6, 9.5, 1.2 };
-  Float_t rangesMax[] = { 44, 6.2, 12, 1.9 };
+  Float_t rangesMin[] = { 136.1, 17.45, 18, 2.4 };
+  Float_t rangesMax[] = { 139.4, 18.95, 23.5, 3.6 };
   
-  if (twoD == 11)
+  for (Int_t i=0; i<2; i++)
   {
-    rangesMin[0] += 1;
-    rangesMax[0] += 1;
+    c->cd(i+1);
+
+    gPad->SetRightMargin(0.05);
+    gPad->SetLeftMargin(0.20);
+    gPad->SetTopMargin(0.15);
+    gPad->SetBottomMargin(0.15);
+
+    if (i == 0)
+    {
+      gPad->SetBottomMargin(0.01);
+    }
+    else if (i == 1)
+    {
+      gPad->SetTopMargin(0.01);
+      gPad->SetBottomMargin(0.7);
+    }
   }
+
+  Int_t i = 0;
+
+  TF1* flowFunc = 0;
   
-  for (Int_t i=0; i<4; i++)
+  c->cd(i+1);
+//   hist = MACHCone(fileName, i, "", kTRUE, 0, &flowFunc);
+  hist = (TH1*) MACHCone(fileName, i, 0, kFALSE, 0)->Clone("hist");
+  hist->GetYaxis()->SetRangeUser(rangesMin[i], rangesMax[i]);
+  hist->GetYaxis()->SetTitleOffset(1.7);
+  hist->Draw();
+//     continue;
+  
+  c->cd(i+1);
+  //c->cd(i+2);
+//   hist2 = MACHCone(fileName, i, 0, kFALSE, 11);
+  hist2 = MACHCone(fileName, i, "SAME", kTRUE, 11, &flowFunc);
+  //hist2->GetYaxis()->SetRangeUser(rangesMin[i], rangesMax[i]);
+//   hist2->Scale(2);
+  hist2->SetLineColor(4);
+  //hist2->Draw("SAME");
+//    return;
+
+  // sample flowFunc with same binning
+  flowFuncHist = (TH1*) hist->Clone("flowFuncHist");
+  flowFuncHist->Reset();
+  
+  for (Int_t i=1; i<=flowFuncHist->GetNbinsX(); i++)
+  {
+    flowFuncHist->SetBinContent(i, flowFunc->Integral(flowFuncHist->GetXaxis()->GetBinLowEdge(i), flowFuncHist->GetXaxis()->GetBinUpEdge(i)));
+    flowFuncHist->SetBinError(i, 0);
+  }
+  
+  i = 0;
+  
+  flowFuncHist->Scale(1.0 / flowFuncHist->GetBinWidth(1));
+  
+//   new TCanvas;
+//   flowFuncHist->Draw();
+//   flowFunc->Draw("SAME");
+//   return;
+
+  // residuals
+  residuals = (TH1*) hist->Clone("res1");
+  residuals->SetYTitle("Residuals");
+  residuals->Divide(flowFuncHist);
+  c->cd(i+2);
+  residuals->Draw();
+  residuals->GetYaxis()->SetNdivisions(505);
+  residuals->GetYaxis()->SetRangeUser(0.995, 1.005);
+
+  residuals = (TH1*) hist2->Clone("res2");
+  residuals->Divide(flowFuncHist);
+  c->cd(i+2);
+  residuals->SetLineColor(4);
+  residuals->Draw("SAME");
+  
+  gPad->SetGridy();
+
+  if (i == 0)
   {
+    legend = new TLegend(0.43, 0.66, 0.91, 0.82);
+    legend->SetFillColor(0);
+    legend->SetTextSize(0.034);
+    legend->SetTextAlign(22);
+    
+    legend->SetHeader("Centrality 0-1%");
+    legend->AddEntry(hist, "|#eta| < 0.8 & All #Delta#eta", "L");
+    legend->AddEntry(hist2, "|#eta| < 0.8 & |#Delta#eta| > 0.8 (2x)", "L");
+    
     c->cd(i+1);
-    MACHCone(fileName, i, "", kTRUE, twoD)->GetYaxis()->SetRangeUser(rangesMin[i], rangesMax[i]);
+    legend->Draw();
   }
   
-  c->SaveAs(Form("machcone_%d.png", twoD));
-  c->SaveAs(Form("machcone_%d.eps", twoD));
-  c->SaveAs(Form("machcone_%d.C", twoD));
+  c->SaveAs(Form("machcone_%d.png", 0));
+  c->SaveAs(Form("machcone_%d.eps", 0));
+  c->SaveAs(Form("machcone_%d.C", 0));
   
   return;
   
@@ -7135,41 +7574,96 @@ void PlotQA(const char* fileName)
   
   TFile::Open(fileName);
   
+  // phys sel
+  Int_t runNumber = 0;
+  list = (TList*) gFile->Get("PhysSel");
+  if (list)
+  {
+    physSel = (AliPhysicsSelection*) list->FindObject("AliPhysicsSelection");
+//     runNumber = physSel->GetCurrentRun();
+  }
+  
+  TString tmp(fileName);
+  tmp.ReplaceAll("alien:///alice/cern.ch/user/j/jgrosseo/gridjob/dir_", "");
+  tmp.ReplaceAll(".root", ".png");
+  tmp.ReplaceAll("/", "-");
+  TString title;
+  title.Form("QA_%d_%s", runNumber, tmp.Data());
+  c = new TCanvas(title, title, 1200, 800);
+  c->Divide(3, 3);
+
   // QA task
   list = (TList*) gFile->Get("histosPhiCorrelationsQA");
+  if (list)
+  {
+    prim = (TH2*) list->FindObject("fDCAPrimaries");
+    dcaxy = prim->ProjectionX("dcaxy", prim->GetYaxis()->FindBin(-3.2), prim->GetYaxis()->FindBin(3.2));
+    dcaz = prim->ProjectionY("dcaz", prim->GetXaxis()->FindBin(-2.4), prim->GetXaxis()->FindBin(2.4));
+    centrCorr = (TH2*) list->FindObject("fCentralityCorrelation");
+    
+    c->cd(1); dcaxy->Draw(); dcaz->SetLineColor(2); dcaz->Draw("SAME");  gPad->SetLogy(); 
+    c->cd(6); centrCorr->Draw("COLZ"); gPad->SetLogz();
   
-  prim = (TH2*) list->FindObject("fDCAPrimaries");
-  dcaxy = prim->ProjectionX("dcaxy", prim->GetYaxis()->FindBin(-3.2), prim->GetYaxis()->FindBin(3.2));
-  dcaz = prim->ProjectionY("dcaz", prim->GetXaxis()->FindBin(-2.4), prim->GetXaxis()->FindBin(2.4));
-  centrCorr = (TH2*) list->FindObject("fCentralityCorrelation");
-  
-  cuts = (AliESDtrackCuts*) list->FindObject("cuts_quality_dca");
-  cluster = cuts->GetNClustersTPC(1);
-  
-  ptall = (TH1F*) cuts->GetPtHist(1)->Clone("ptall");
-  if (ptall->Integral() > 0)
-    ptall->Scale(1.0 / ptall->Integral());
-  
-  check_its = (AliESDtrackCuts*) list->FindObject("check_its");
-  ptIts = (TH1F*) check_its->GetPtHist(1)->Clone("ptIts");
-  if (ptIts->Integral() > 0)
-    ptIts->Scale(1.0 / ptIts->Integral());
-  
-  global_cuts = (AliESDtrackCuts*) list->FindObject("global_cuts");
-  ptItsAcc = (TH1*) global_cuts->GetPtHist(1)->Clone("ptItsAcc");
-  if (ptItsAcc->Integral() > 0)
-    ptItsAcc->Scale(1.0 / ptItsAcc->Integral());
+    cuts = (AliESDtrackCuts*) list->FindObject("cuts_quality_dca");
+    if (cuts)
+    {
+      cluster = cuts->GetNClustersTPC(1);
+      c->cd(3); cluster->Draw();
+    
+      ptall = (TH1F*) cuts->GetPtHist(1)->Clone("ptall");
+      if (ptall->Integral() > 0)
+       ptall->Scale(1.0 / ptall->Integral());
+    
+      c->cd(7); 
+      ptall->Draw(); 
+
+      TH1* ptIts = 0;
+      check_its = (AliESDtrackCuts*) list->FindObject("check_its");
+      if (check_its)
+      {
+       ptIts = (TH1F*) check_its->GetPtHist(1)->Clone("ptIts");
+       if (ptIts->Integral() > 0)
+         ptIts->Scale(1.0 / ptIts->Integral());
+      }
+      
+      TH1* ptItsAcc = 0;
+      global_cuts = (AliESDtrackCuts*) list->FindObject("global_cuts");
+      if (global_cuts)
+      {
+       ptItsAcc = (TH1*) global_cuts->GetPtHist(1)->Clone("ptItsAcc");
+       if (ptItsAcc->Integral() > 0)
+         ptItsAcc->Scale(1.0 / ptItsAcc->Integral());
+      }
+    
+      if (ptIts)
+      {
+       ptIts->SetLineColor(2); 
+       ptIts->Draw("SAME");
+      }
+      
+      if (ptItsAcc)
+      {
+       ptItsAcc->SetLineColor(4); 
+       ptItsAcc->Draw("SAME"); 
+      }
+      
+      gPad->SetLogy();
+    }
+  }
   
   // centrality task
   list = (TList*) gFile->Get("CentralityStat");
-  centrQuality = (TH1*) list->FindObject("fHOutQuality");
-  
-  // phys sel
-  list = (TList*) gFile->Get("PhysSel");
-  physSel = (AliPhysicsSelection*) list->FindObject("AliPhysicsSelection");
+  TH1* centrQuality = 0;
+  if (list)
+  {
+    centrQuality = (TH1*) list->FindObject("fHOutQuality");
+    
+    c->cd(4); if (centrQuality) centrQuality->Draw();
+  }
   
   // phi corr task
   AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+  ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0))->FillParent();
   centr = h->GetCentralityDistribution();
   NormalizeToBinWidth(centr);
   Int_t events = (Int_t) h->GetEventCount()->ProjectionX()->GetBinContent(3);
@@ -7178,30 +7672,32 @@ void PlotQA(const char* fileName)
   dphi_corr = h->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed, 0, 4.01, 14.99, 1, 8);
   
   AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);
+  ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->FillParent();
   hMixed->SetPtRange(1.01, 3.99);
   dphi_corr_mixed = hMixed->GetUEHist(2)->GetUEHist(AliUEHist::kCFStepReconstructed, 0, 4.01, 14.99, 1, 8);
   
-  TString tmp(fileName);
-  tmp.ReplaceAll("alien:///alice/cern.ch/user/j/jgrosseo/gridjob/dir_", "");
-  tmp.ReplaceAll(".root", ".png");
-  tmp.ReplaceAll("/", "-");
-  TString title;
-  title.Form("QA_%d_%s", physSel->GetCurrentRun(), tmp.Data());
-  c = new TCanvas(title, title, 1200, 800);
-  c->Divide(3, 3);
+  if (runNumber != 0 && runNumber != h->GetRunNumber())
+    AliFatal("Run numbers inconsistent");
+
+  Printf("%d", h->GetRunNumber());
+  if (runNumber == 0)
+    runNumber = h->GetRunNumber();
+  c->cd(2); dphi_corr->Draw(); dphi_corr->GetYaxis()->SetRangeUser(dphi_corr->GetMinimum() * 0.9, dphi_corr->GetMaximum() * 1.1); dphi_corr_mixed->SetLineColor(2); dphi_corr_mixed->Draw("SAME");
+  c->cd(5); centr->Draw("HIST");
   
+  c->cd(1);
   latex = new TLatex(0.15, 0.8, Form("%d events", events));
   latex->SetTextSize(0.075);
   latex->SetNDC();
+  latex->Draw();
   
-  c->cd(1); dcaxy->Draw(); dcaz->SetLineColor(2); dcaz->Draw("SAME");  gPad->SetLogy(); latex->Draw();
-  c->cd(2); dphi_corr->Draw(); dphi_corr->GetYaxis()->SetRangeUser(dphi_corr->GetMinimum() * 0.9, dphi_corr->GetMaximum() * 1.1); dphi_corr_mixed->SetLineColor(2); dphi_corr_mixed->Draw("SAME");
-  c->cd(3); cluster->Draw();
-  c->cd(4); centrQuality->Draw();
-  c->cd(5); centr->Draw("HIST");
-  c->cd(6); centrCorr->Draw("COLZ"); gPad->SetLogz();
-  c->cd(7); ptall->Draw(); ptIts->SetLineColor(2); ptIts->Draw("SAME"); ptItsAcc->SetLineColor(4); ptItsAcc->Draw("SAME"); gPad->SetLogy();
+  c->cd(8);
+  h->GetEventCount()->Draw("TEXT");
   
+  title.Form("QA_%d_%s", runNumber, tmp.Data());
+  c->SetTitle(title);
+
   c->SaveAs(Form("qa/%s", c->GetTitle()));
 }
 
@@ -7604,3 +8100,127 @@ void NormalizeTo(TGraphErrors* graph, Float_t normalizeTo)
        }       
 }
 
+void CompareMixedEvent(const char* fileName)
+{
+  loadlibs();
+  
+  Float_t leadingPtArr[] = { 6.0, 8.0, 10.0, 15.0, 15.0 };
+  Float_t assocPtArr[] =     { 0.5, 1.5, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0 };
+  Int_t leadingPtOffset = 2;
+  Int_t centralityBins[] = { 0, 0, 1, 6, 9, 16 };
+
+  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);  
+  
+  Int_t i = 1;
+  for (Int_t j=2; j<5; j++)
+  {
+    gpTMin = assocPtArr[j] + 0.01;
+    gpTMax = assocPtArr[j+1] - 0.01;
+    
+    gpTMin = 3.0;
+    gpTMax = 6.0;
+
+    SetupRanges(hMixed);
+    
+    TH2* mixed = hMixed->GetUEHist(2)->GetUEHist(6, 0, leadingPtArr[i] + 0.01, leadingPtArr[i+leadingPtOffset] - 0.01, centralityBins[j], centralityBins[j+1]-1, 1);
+  
+    // compare deta
+    
+    TH1* histMixedproj = mixed->ProjectionY();
+    histMixedproj->Scale(1.0 / mixed->GetNbinsX());
+    
+    for (Int_t x=1; x<=mixed->GetNbinsX(); x++)
+      for (Int_t y=1; y<=mixed->GetNbinsY(); y++)
+       mixed->SetBinContent(x, y, histMixedproj->GetBinContent(y));
+
+    histMixedproj->Scale(1.0 / (0.5 * (histMixedproj->GetBinContent(histMixedproj->GetXaxis()->FindBin(-0.01)) + histMixedproj->GetBinContent(histMixedproj->GetXaxis()->FindBin(0.01)))));
+      
+    histMixedproj->DrawCopy((j == 2) ? "" : "SAME")->SetLineColor(j-1);
+  }
+}
+
+void FillParentTHnSparse(const char* fileName)
+{
+  loadlibs();
+
+  TList* list = 0;
+  
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName, &list);
+  ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0))->FillParent();
+  ((AliTHn*) h->GetUEHist(2)->GetTrackHist(0))->DeleteContainers();
+  
+  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);  
+  ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->FillParent();
+  ((AliTHn*) hMixed->GetUEHist(2)->GetTrackHist(0))->DeleteContainers();
+  
+  TString newFileName(fileName);
+  newFileName.ReplaceAll(".root", "");
+  newFileName += "_.root";
+
+  file3 = TFile::Open(newFileName, "RECREATE");
+  file3->mkdir("PWG4_PhiCorrelations");
+  file3->cd("PWG4_PhiCorrelations");
+  list->Write("histosPhiCorrelations", TObject::kSingleKey);
+  file3->Close();
+}
+
+void CompareZVertex(const char* fileName)
+{
+  loadlibs();
+  
+  AliUEHistograms* h = (AliUEHistograms*) GetUEHistogram(fileName);
+  AliUEHistograms* hMixed = (AliUEHistograms*) GetUEHistogram(fileName, 0, kTRUE);  
+  
+  axis = h->GetUEHist(2)->GetEventHist()->GetAxis(2, 0);
+  
+  gpTMin = 0.151;
+  gpTMax = 1.0;
+  
+  SetupRanges(h);
+  SetupRanges(hMixed);
+  
+  TFile::Open("test.root", "RECREATE");
+       
+  for (Int_t i=0; i<=axis->GetNbins(); i++)
+  {
+    TH1* hist = 0;
+    if (i > 0)
+    {
+      h->SetZVtxRange(h->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetBinLowEdge(i) + 0.01, h->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetBinUpEdge(i) - 0.01);
+      hMixed->SetZVtxRange(h->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetBinLowEdge(i) + 0.01, h->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetBinUpEdge(i) - 0.01);
+    }
+      
+    GetDistAndFlow(h, hMixed, &hist, 0, 6, 0, 10, 0.51, 1.99, 1, kTRUE, 0, kFALSE);
+    
+    new TCanvas;
+    hist->DrawCopy("SURF1");
+    
+    hist->Write(Form("detadphi_%d", i));
+    
+//     if (i == 0)   break;
+
+    continue;
+    
+    hist->SetLineColor(i+1);
+    hist->Scale(1.0 / hist->Integral() / hist->GetBinWidth(1));
+    hist->Draw((i == 0) ? "" : "SAME");
+  }
+  
+  gFile->Close();
+}
+
+void DrawZRanges(Float_t min, Float_t max)
+{
+  TFile::Open("test.root");
+  
+  for (Int_t i=0; i<8; i++)
+  {
+    hist = (TH2*) gFile->Get(Form("detadphi_%d", i));
+    
+    proj = hist->ProjectionY("proj", hist->GetXaxis()->FindBin(min), hist->GetXaxis()->FindBin(max));
+    proj->Scale(1.0 / (hist->GetXaxis()->FindBin(max) - hist->GetXaxis()->FindBin(min) + 1));
+    
+    proj->SetLineColor(i+1);
+    proj->DrawCopy((i == 0) ? "" : "SAME");
+  }
+}