From 7e70f6e0b0d5c9ec332b1ebf6fe405ab55921d5a Mon Sep 17 00:00:00 2001 From: pcrochet Date: Sun, 22 Jul 2012 13:53:23 +0000 Subject: [PATCH] Add task for single muon flow estimation with the EP method (Diego) --- PWG/CMakelibPWGmuon.pkg | 1 + PWG/PWGmuonLinkDef.h | 1 + PWG/muon/AddTaskFlowSingleMu.C | 51 ++ PWG/muon/AliAnalysisTaskFlowSingleMu.cxx | 796 +++++++++++++++++++++++ PWG/muon/AliAnalysisTaskFlowSingleMu.h | 80 +++ 5 files changed, 929 insertions(+) create mode 100644 PWG/muon/AddTaskFlowSingleMu.C create mode 100644 PWG/muon/AliAnalysisTaskFlowSingleMu.cxx create mode 100644 PWG/muon/AliAnalysisTaskFlowSingleMu.h diff --git a/PWG/CMakelibPWGmuon.pkg b/PWG/CMakelibPWGmuon.pkg index c9f50b74b97..73e2c79a975 100644 --- a/PWG/CMakelibPWGmuon.pkg +++ b/PWG/CMakelibPWGmuon.pkg @@ -61,6 +61,7 @@ set ( SRCS muon/AliMuonPairCuts.cxx muon/AliMergeableCollection.cxx muon/AliVAnalysisMuon.cxx + muon/AliAnalysisTaskFlowSingleMu.cxx ) string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) diff --git a/PWG/PWGmuonLinkDef.h b/PWG/PWGmuonLinkDef.h index 628c3d457d5..66e18c15ae2 100644 --- a/PWG/PWGmuonLinkDef.h +++ b/PWG/PWGmuonLinkDef.h @@ -42,5 +42,6 @@ #pragma link C++ class AliMergeableCollection+; #pragma link C++ class AliMergeableCollectionIterator+; #pragma link C++ class AliVAnalysisMuon+; +#pragma link C++ class AliAnalysisTaskFlowSingleMu+; #endif diff --git a/PWG/muon/AddTaskFlowSingleMu.C b/PWG/muon/AddTaskFlowSingleMu.C new file mode 100644 index 00000000000..ab95a5ea787 --- /dev/null +++ b/PWG/muon/AddTaskFlowSingleMu.C @@ -0,0 +1,51 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) +#include "TString.h" +#include "TObjArray.h" + +#include "AliLog.h" +#include "AliVEventHandler.h" + +#include "AliAnalysisManager.h" +#include "AliAnalysisDataContainer.h" + +#include "AliMuonTrackCuts.h" +#include "AliAnalysisTaskFlowSingleMu.h" +#endif + +AliAnalysisTaskFlowSingleMu* AddTaskFlowSingleMu(Bool_t isMC = kFALSE) +{ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddtaskFlowSingleMu", "No analysis manager to connect to."); + return NULL; + } + + TString type = mgr->GetInputEventHandler()->GetDataType(); + if (!type.Contains("ESD") && !type.Contains("AOD")) { + ::Error("AddtaskFlowSingleMu", "FlowSingleMu task needs the manager to have an ESD or AOD input handler."); + return NULL; + } + + // Create container + TString currName = ""; + TString outputfile = mgr->GetCommonFileName(); + if ( ! outputfile.IsNull() ) outputfile += ":PWGHF_FlowSingleMu"; + else outputfile = "FlowSingleMuAnalysis.root"; + + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("FlowSingleMuOut",TObjArray::Class(),AliAnalysisManager::kOutputContainer,outputfile); + + // Create cuts + AliMuonTrackCuts* muonTrackCuts = new AliMuonTrackCuts("TestStandardMuonTrackCuts", "TestStandardMuonTrackCuts"); + muonTrackCuts->SetIsMC(isMC); + + // Create task + AliAnalysisTaskFlowSingleMu *flowSingleMuTask = new AliAnalysisTaskFlowSingleMu("FlowSingleMuTask", *muonTrackCuts); + if ( isMC ) flowSingleMuTask->SetTrigClassPatterns("ANY"); + mgr->AddTask(flowSingleMuTask); + + // Connect containers + mgr->ConnectInput (flowSingleMuTask, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput (flowSingleMuTask, 1, coutput1); + + return flowSingleMuTask; +} diff --git a/PWG/muon/AliAnalysisTaskFlowSingleMu.cxx b/PWG/muon/AliAnalysisTaskFlowSingleMu.cxx new file mode 100644 index 00000000000..7a9bced7a75 --- /dev/null +++ b/PWG/muon/AliAnalysisTaskFlowSingleMu.cxx @@ -0,0 +1,796 @@ +/************************************************************************** + * Copyright(c) 1998-2007, 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: AliAnalysisTaskFlowSingleMu.cxx 55545 2012-04-04 07:16:39Z pcrochet $ */ + +//----------------------------------------------------------------------------- +/// \class AliAnalysisTaskFlowSingleMu +/// Analysis task for flow of single muons in the spectrometer. +/// The output is a list of histograms and CF containers. +/// The macro class can run on AODs or ESDs. +/// If Monte Carlo information is present, some basics checks are performed. +/// +/// \author Diego Stocco +//----------------------------------------------------------------------------- + +#define AliAnalysisTaskFlowSingleMu_cxx + +#include "AliAnalysisTaskFlowSingleMu.h" + +// ROOT includes +#include "TROOT.h" +#include "TH1.h" +#include "TH2.h" +#include "TGraphErrors.h" +#include "TAxis.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TMath.h" +#include "TObjString.h" +#include "TObjArray.h" +#include "TF1.h" +#include "TStyle.h" +//#include "TMCProcess.h" +#include "TFitResultPtr.h" +#include "TFile.h" +#include "TArrayI.h" +#include "TArrayD.h" +#include "TSystem.h" +#include "TGraphErrors.h" +#include "TLatex.h" + +// STEER includes +#include "AliAODEvent.h" +#include "AliAODTrack.h" +#include "AliAODMCParticle.h" +#include "AliMCEvent.h" +#include "AliMCParticle.h" +#include "AliESDEvent.h" +#include "AliESDMuonTrack.h" +#include "AliVHeader.h" +#include "AliAODMCHeader.h" +#include "AliStack.h" +#include "AliEventplane.h" + +// ANALYSIS includes +#include "AliAnalysisManager.h" + +// CORRFW includes +#include "AliCFContainer.h" +#include "AliCFGridSparse.h" +#include "AliCFEffGrid.h" + +// PWG includes +#include "AliVAnalysisMuon.h" +#include "AliMergeableCollection.h" +#include "AliCounterCollection.h" +#include "AliMuonTrackCuts.h" + + +/// \cond CLASSIMP +ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context +/// \endcond + + +//________________________________________________________________________ +AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() : + AliVAnalysisMuon(), + fEPKeys(0x0), + fRandom(0x0) +{ + /// Default ctor. +} + +//________________________________________________________________________ +AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) : + AliVAnalysisMuon(name, cuts), + fEPKeys(0x0), + fRandom(0x0) +{ + // + /// Constructor. + // + TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4"; + fEPKeys = epMethods.Tokenize(" "); +} + + +//________________________________________________________________________ +AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu() +{ + // + /// Destructor + // + delete fEPKeys; + delete fRandom; +} + +//___________________________________________________________________________ +void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects() +{ + // + /// Create outputs + // + + if ( ! fRandom ) fRandom = new TRandom3(); + // Set seed here, to avoid having the same seed for all jobs in grid + fRandom->SetSeed(0); + + Int_t nPtBins = 160; + Double_t ptMin = 0., ptMax = 80.; + TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c"); + + Int_t nEtaBins = 25; + Double_t etaMin = -4.5, etaMax = -2.; + TString etaName("Eta"), etaTitle("#eta"), etaUnits(""); + + Int_t nPhiBins = 36; + Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi(); + TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad"); + + Int_t nDeltaPhiBins = 18; + Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi(); + TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad"); + + Int_t nChargeBins = 2; + Double_t chargeMin = -2, chargeMax = 2.; + TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e"); + + Int_t nMotherTypeBins = kNtrackSources; + Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5; + TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits(""); + + Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins}; + Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin}; + Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax}; + TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle}; + TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits}; + + Int_t iobj=0; + TH1* histo = 0x0; + TString currTitle = ""; + for ( Int_t iep=0; iepGetEntries(); iep++ ) { + TString currEP = fEPKeys->At(iep)->GetName(); + AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins); + + TString currName = ""; + for ( Int_t idim = 0; idimSetVarTitle(idim, currName.Data()); + cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]); + } + + TString stepTitle[kNsteps] = {"reconstructed", "generated"}; + + TAxis* currAxis = 0x0; + for (Int_t istep=0; istepSetStepTitle(istep, stepTitle[istep].Data()); + AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep); + + currAxis = gridSparse->GetAxis(kHvarMotherType); + for ( Int_t ibin=0; ibinGetEntries(); ibin++ ) { + currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName()); + } + } + + AddObjectToCollection(cfContainer, iobj++); + + Bool_t isTPC = currEP.Contains("TPC"); + Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.; + Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.; + + histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane); + histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data())); + AddObjectToCollection(histo, iobj++); + + for ( Int_t jep=iep+1; jepGetEntries(); jep++ ) { + TString auxEP = fEPKeys->At(jep)->GetName(); + currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data()); + histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.); + histo->GetXaxis()->SetTitle(currTitle.Data()); + AddObjectToCollection(histo, iobj++); + } + + currTitle = Form("cos(#psi_{1} - #psi_{2})"); + histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.); + histo->GetXaxis()->SetTitle(currTitle.Data()); + + AddObjectToCollection(histo, iobj++); + } // loop on EPs + + fMuonTrackCuts->Print("mask"); +} + +//________________________________________________________________________ +void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality) +{ + // + /// Fill output objects + // + + // + // Base event cuts + // + if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return; + //Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes + //if ( isPileupFromSPD ) return; + if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) return; + + TArrayD psiPlane(fEPKeys->GetEntriesFast()); + psiPlane.Reset(-999.); + const Double_t kLowestVal = -10.; + Int_t harmonic = 2; + + // + // Fill EP info + // + AliEventplane* eventPlane = InputEvent()->GetEventplane(); + psiPlane[0] = eventPlane->GetEventplane("V0A",InputEvent(),harmonic); + psiPlane[1] = eventPlane->GetEventplane("V0C",InputEvent(),harmonic); + psiPlane[2] = eventPlane->GetEventplane("Q"); + if ( psiPlane[2] < 0.) psiPlane[2] = -999.; + Double_t qx = 0., qy = 0.; + psiPlane[3] = eventPlane->CalculateVZEROEventPlane(InputEvent(),0,harmonic,qx,qy); // V0C ring 1 + psiPlane[4] = eventPlane->CalculateVZEROEventPlane(InputEvent(),3,harmonic,qx,qy); // V0C ring 3 + //psiPlane[3] = fRandom->Uniform(-TMath::Pi()/2., TMath::Pi()/2.); + + Bool_t noValidEP = kTRUE; + TArrayD psiPlaneCalc(fEPKeys->GetEntriesFast()); + for ( Int_t iep=0; iepGetEntriesFast(); iep++ ) { + psiPlaneCalc[iep] = psiPlane[iep]; + if ( psiPlane[iep] < kLowestVal ) continue; + if ( psiPlane[iep] < 0. ) psiPlaneCalc[iep] += TMath::Pi(); + noValidEP = kFALSE; + } + if ( noValidEP ) return; + + for ( Int_t iep=0; iepGetEntriesFast(); iep++ ) { + if ( psiPlane[iep] < kLowestVal ) continue; + for ( Int_t itrig=0; itrigGetString(); + TString currEP = fEPKeys->At(iep)->GetName(); + ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hEventPlane%s",currEP.Data())))->Fill(psiPlane[iep]); + for ( Int_t jep=iep+1; jepGetEntriesFast(); jep++ ) { + if ( psiPlane[jep] < kLowestVal ) continue; + ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub3%s%s",currEP.Data(), fEPKeys->At(jep)->GetName())))->Fill(TMath::Cos(2.*(psiPlaneCalc[iep]-psiPlaneCalc[jep]))); + } // loop on auxialiary EP + if ( iep == 2 ) { + TVector2* qsub1 = eventPlane->GetQsub1(); + TVector2* qsub2 = eventPlane->GetQsub2(); + if ( qsub1 && qsub2 ) + ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.)); + } + } // loop on trigger + } // loop on EP + + // + // Fill track info + // + Double_t containerInput[kNvars]; + AliVParticle* track = 0x0; + for ( Int_t istep = 0; istep<2; ++istep ) { + Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks(); + for (Int_t itrack = 0; itrack < nTracks; itrack++) { + track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack); + + Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 ); + if ( ! isSelected ) continue; + + Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track); + + + containerInput[kHvarPt] = track->Pt(); + containerInput[kHvarEta] = track->Eta(); + containerInput[kHvarPhi] = track->Phi(); + containerInput[kHvarCharge] = track->Charge()/3.; + containerInput[kHvarMotherType] = (Double_t)trackSrc; + + for ( Int_t iep=0; iepGetEntriesFast(); iep++ ) { + if ( psiPlane[iep] < kLowestVal ) continue; + TString currEP = fEPKeys->At(iep)->GetName(); + Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep]; + + // The difference between phi - psiPlane must be between (0, pi) + // Since the reaction plain is symmetric in pi, + // we can do in such a way that: + // -pi < delta_phi + N*pi < pi + // We choose N accrodingly + Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() ); + deltaPhi += nPi * TMath::Pi(); + containerInput[kHvarDeltaPhi] = deltaPhi; + + for ( Int_t itrig=0; itrigGetString(); + if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue; + ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep); + } // loop on selected trigger classes + } // loop on event plane + } // loop on tracks + } // loop on container steps +} + + +//________________________________________________________________________ +TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange) +{ + // + /// Get the range from string + // + TArrayD centralityRange(2); + centralityRange.Reset(-999.); + TObjArray* array = sRange.Tokenize("_"); + if ( array->GetEntries() == 2 ) { + for ( Int_t ientry=0; ientry<2; ientry++ ) { + centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof(); + } + } + delete array; + return centralityRange; +} + +//________________________________________________________________________ +void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) { + // + /// Draw some histograms at the end. + // + + AliVAnalysisMuon::Terminate(""); + + if ( ! fMergeableCollection ) return; + + TString physSel = fTerminateOptions->At(0)->GetName(); + TString trigClassName = fTerminateOptions->At(1)->GetName(); + TString centralityRange = fTerminateOptions->At(2)->GetName(); + TString furtherOpt = fTerminateOptions->At(3)->GetName(); + + TObjArray* centralityClasses = centralityRange.Tokenize(","); + TArrayD centralityArray(centralityClasses->GetEntries()+1); + Int_t ib = 0; + for ( Int_t icent=0; icentGetEntries(); icent++ ) { + TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName()); + if ( icent == 0 ) centralityArray.AddAt(range[0],ib++); + centralityArray.AddAt(range[1],ib++); + } + + TString selectEP = ""; + TObjArray resoEParray; + resoEParray.SetOwner(); + + TString outFilename = ""; + TObjArray* optArr = furtherOpt.Tokenize(" "); + TString currName = ""; + TString resoFilename = ""; + for ( Int_t iopt=0; ioptGetEntries(); iopt++ ) { + currName = optArr->At(iopt)->GetName(); + if ( currName.Contains(".root") ) outFilename = currName; + else if ( currName.Contains(".txt") ) resoFilename = currName; + else { + for ( Int_t iep=0; iepGetEntries(); iep++ ) { + if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data())); + } + } + } + delete optArr; + + if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName(); + + furtherOpt.ToUpper(); + + Bool_t v2only = furtherOpt.Contains("V2ONLY"); + + TAxis* customBins = 0x0; + + //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.}; + Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE + //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.}; + //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.}; + //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.}; + Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1; + customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE + + TCanvas *can = 0x0, *showCan = 0x0; + Int_t xshift = 100; + Int_t yshift = 20; + Int_t igroup1 = 0; + Int_t igroup2 = 0; + + // + /// Read plane resolution file (if present) + // + TObjArray resoCentrality; + resoCentrality.SetOwner(); + TArrayD resolution[3]; + for ( Int_t ires=0; ires<3; ires++ ) { + resolution[ires].Set(50); + } + + Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) ); + Bool_t internalReso = kTRUE; + if ( externalReso ) { + // If plane resolution file exists, fill resolution + TString currLine = ""; + ifstream inFile(resoFilename.Data()); + if (inFile.is_open()) { + ib = 0; + while (! inFile.eof() ) { + currLine.ReadLine(inFile,kFALSE); + if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue; + TObjArray* array = currLine.Tokenize(" "); + resoCentrality.Add(new TObjString(array->At(0)->GetName())); + for ( Int_t ires=0; ires<3; ires++ ) { + resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib); + } + delete array; + ib++; + } // ! EOF + } // file is open + for ( Int_t ires=0; ires<3; ires++ ) { + resolution[ires].Set(ib); + } + } // file exists + else { + // If no plane resolution file exist, + // use the centralities in terminateOption + // and set the resolution to 1 and errors to 0. + + for ( Int_t ires=0; ires<3; ires++ ) { + resolution[ires].Set(centralityClasses->GetEntries()); + Double_t resetVal = ( ires == 0 ) ? 1. : 0.; + resolution[ires].Reset(resetVal); + } + TArrayD currCos(3); + for ( Int_t icent=0; icentGetEntries(); icent++ ) { + resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName())); + Int_t icos = 0; + Double_t sumRelErrSquare = 0.; + for ( Int_t isub=0; isub(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName()))); + if ( ! histo || histo->Integral() == 0 ) continue; + currCos[icos] = histo->GetMean(); + if ( currCos[icos] <= 0. ) continue; + Double_t relErr = histo->GetMeanError() / currCos[icos]; + sumRelErrSquare += relErr * relErr; + icos++; + } // loop on sub events + } // loop on sub events + if ( icos != 3 ) { + printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName()); + internalReso = kFALSE; + continue; + } + resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2])); + resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare); + printf("Resolution %s %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]); + } + } + + Bool_t hasResolution = ( externalReso || internalReso ); + + + TArrayD resoCentrArray(resoCentrality.GetEntries()+1); + ib = 0; + for ( Int_t icent=0; icentGetName()); + if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++); + resoCentrArray.AddAt(range[1], ib++); + } + + if ( hasResolution ) { + currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName()); + can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); + TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries()); + for ( Int_t icent=0; icentSetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]); + resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]); + } + resoGraph->SetTitle(Form("Resolution %s", selectEP.Data())); + resoGraph->GetXaxis()->SetTitle("Centrality"); + resoGraph->GetYaxis()->SetTitle("Resolution"); + resoGraph->Draw("apz"); + } + + const Int_t kNresoCentr = resoCentrality.GetEntries(); + + // + // Create gridSparse array + // + AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr]; + TString stepName[2]; + for ( Int_t icent=0; icent ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) ); + if ( ! cfContainer ) continue; + for ( Int_t istep=0; istepGetStepTitle(istep); + gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep); + } // loop on step + } // loop on centrality + + Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange}; + + Bool_t isMC = furtherOpt.Contains("MC"); + Int_t firstSrc = ( isMC ) ? 0 : kUnidentified; + Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified; + if ( ! isMC ) srcColors[kUnidentified] = 1; + + TList outList; + + TString graphTypeName[3] = {"raw","correctStat","correctSyst"}; + Int_t drawOrder[3] = {0, 2, 1}; + TString drawOpt[3] = {"apz","pz","a2"}; + Int_t nTypes = ( hasResolution ) ? 3 : 1; + + TString histoName = "", histoPattern = ""; + /////////// + // Flow // + /////////// + gStyle->SetOptFit(1111); + TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))"; + // + TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi()); + if ( v2only ) func->SetParNames("scale", "v2"); + else func->SetParNames("scale","v1", "v2", "v3"); + Int_t v2par = ( v2only ) ? 1 : 2; + + for ( Int_t istep=0; istepAt(isrc)->GetName()); + TGraphErrors* flowVsCentrality[3]; + for ( Int_t itype=0; itypeSetName(histoName.Data()); + flowVsCentrality[itype]->SetTitle(histoName.Data()); + } + for ( Int_t icent=0; icentGetEntries(); icent++ ) { + TH2* hDeltaPhi2D = 0x0; + TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName()); + TArrayD weightedReso(3); + weightedReso.Reset(0.); + Double_t nMuons = 0.; + TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName()); + for ( Int_t rcent=0; rcent centralityArray[icent+1] ) continue; + AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent]; + if ( ! gridSparse || gridSparse->GetEntries() == 0. ) continue; + + SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501); + SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN"); + TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt); + + Double_t currYield = auxHisto->Integral(); + if ( currYield == 0. ) continue; + nMuons += currYield; + debugString += Form("\nAdding %s yield %g reso", resoCentrality.At(rcent)->GetName(), currYield); + for ( Int_t ires=0; ires<3; ires++ ) { + weightedReso[ires] += resolution[ires][rcent] * currYield; + debugString += Form(" %g", resolution[ires][rcent]); + } + histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data()); + if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast(auxHisto->Clone(histoName.Data())); + else hDeltaPhi2D->Add(auxHisto); + delete auxHisto; + } + + if ( ! hDeltaPhi2D ) continue; + + debugString += Form("\nWeighted reso "); + for ( Int_t ires=0; ires<3; ires++ ) { + weightedReso[ires] = weightedReso[ires]/nMuons; + debugString += Form(" %g", weightedReso[ires]); + } + AliDebug(1,debugString.Data()); + + TAxis* ptAxis = hDeltaPhi2D->GetYaxis(); + + Int_t ipad = 0; + TGraphErrors* flowVsPt[3]; + for ( Int_t itype=0; itypeSetName(histoName.Data()); + flowVsPt[itype]->SetTitle(histoName.Data()); + } + + if ( ! customBins ) customBins = static_cast(ptAxis->Clone()); + + + for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) { + Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt; + Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt; + Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4); + Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4); + Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin); + Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin); + TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin); + projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax)); + if ( projHisto->Integral() < 50 ) break; + if ( ipad%4 == 0 ) { + currName = projHisto->GetName(); + currName.Append("_can"); + showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); + showCan->Divide(2,2); + ipad = 0; + } + showCan->cd(++ipad); + if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01); + else { + func->SetParameters(projHisto->Integral(), 0., 0.01, 0.); + func->FixParameter(1,0.); + } + projHisto->Fit(func,"R"); + for ( Int_t itype=0; itypeGetParameter(v2par); + Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par); + Double_t yVal = rawVal/resoVal; + Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin); + Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin); + Int_t ipoint = currGraph->GetN(); + currGraph->SetPoint(ipoint,xVal,yVal); + if ( itype > 0 && ipt != 0 ) { + // For the plot vs. pt the resolution error is fully correlated. + // Hence we do not use add it bin by bin, but rather write it + // on the title + resoErr = 0.; + TString graphTitle = currGraph->GetTitle(); + if ( ! graphTitle.Contains("|") ) { + graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]); + currGraph->SetTitle(graphTitle.Data()); + } + } + Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.; + Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.; + Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel); + currGraph->SetPointError(ipoint,xErr,yErr); + } // loop on type + } // loop on pt bins + for ( Int_t itype=0; itypeGetName(); + currName.Append("Can"); + can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); + igroup2++; + can->cd(); + } + flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle()); + flowVsPt[currType]->GetYaxis()->SetTitle("v2"); + flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25); + flowVsPt[currType]->SetFillStyle(0); + flowVsPt[currType]->Draw(drawOpt[currType].Data()); + outList.Add(flowVsPt[currType]); + } // loop on types + } // loop on centrality + for ( Int_t itype=0; itypeGetN() == 0 ) continue; + if ( itype < 2 ) { + currName = flowVsCentrality[currType]->GetName(); + currName.Append("Can"); + can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); + can->cd(); + } + flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality"); + flowVsCentrality[currType]->GetYaxis()->SetTitle("v2"); + flowVsCentrality[currType]->SetFillStyle(0); + flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25); + flowVsCentrality[currType]->Draw(drawOpt[currType].Data()); + outList.Add(flowVsCentrality[currType]); + igroup2++; + } // loop on type + igroup1++; + igroup2 = 0; + } // loop on track sources + } // loop on steps + + delete customBins; + + + /////////////////////// + // Event plane check // + /////////////////////// + igroup1++; + igroup2 = 0; + Int_t icolor=0; + currName ="eventPlaneDistrib"; + can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600); + can->SetGridy(); + igroup2++; + + Bool_t addOffsetToCentrality = kFALSE; + + TLegend* leg = new TLegend(0.7,0.7,0.92,0.92); + leg->SetBorderSize(1); + for ( Int_t icent=0; icentGetEntries(); icent++ ) { + TH1* histo = static_cast ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) ); + if ( ! histo ) continue; + histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName())); + if ( histo->Integral() < 50. ) continue; + if ( histo->GetSumw2N() == 0 ) histo->Sumw2(); + histo->Fit("pol0","R0"); + TF1* pol0 = histo->GetFunction("pol0"); + histo->Scale(1./pol0->GetParameter(0)); + Int_t offset = ( addOffsetToCentrality ) ? icolor : 0; + TGraphErrors* graph = new TGraphErrors(); + graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data())); + for ( Int_t ipoint=0; ipointGetXaxis()->GetNbins(); ipoint++ ) { + Int_t ibin = ipoint+1; + graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset); + graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin)); + } + graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle()); + graph->GetYaxis()->SetTitle("Data/Fit (pol0)"); + TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName()); + if ( addOffsetToCentrality ) { + graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries()); + graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0); + legTitle += Form(" ( + %i)", offset); + } + Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor; + graph->SetLineColor(currColor); + graph->SetMarkerColor(currColor); + graph->SetMarkerStyle(20+icolor); + TString currDrawOpt = "pz"; + if ( icolor == 0 ) currDrawOpt += "a"; + graph->Draw(currDrawOpt.Data()); + legTitle.ReplaceAll("_"," - "); + leg->AddEntry(graph,legTitle.Data(),"lp"); + if ( addOffsetToCentrality ) { + TLatex* latex = new TLatex(histo->GetXaxis()->GetXmax()+0.5*histo->GetXaxis()->GetBinWidth(1), 1.+(Double_t)offset, Form("#chi^{2}/NDF = %.3g", pol0->GetChisquare()/(Double_t)pol0->GetNDF())); + latex->SetTextColor(currColor); + latex->SetTextSize(0.025); + latex->Draw("same"); + } + icolor++; + } + leg->Draw("same"); + + + ////////////////////// + // Event statistics // + ////////////////////// + printf("\nTotal analyzed events:\n"); + TString evtSel = Form("trigger:%s", trigClassName.Data()); + fEventCounters->PrintSum(evtSel.Data()); + printf("Physics selected analyzed events:\n"); + evtSel = Form("trigger:%s/selected:yes", trigClassName.Data()); + fEventCounters->PrintSum(evtSel.Data()); + + TString countPhysSel = "any"; + if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes"; + else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no"; + countPhysSel.Prepend("selected:"); + printf("Analyzed events vs. centrality:\n"); + evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data()); + fEventCounters->Print("centrality",evtSel.Data(),kTRUE); + + if ( ! outFilename.IsNull() ) { + printf("\nWriting output file %s\n", outFilename.Data()); + TFile* file = TFile::Open(outFilename.Data(), "RECREATE"); + outList.Write(); + file->Close(); + } +} diff --git a/PWG/muon/AliAnalysisTaskFlowSingleMu.h b/PWG/muon/AliAnalysisTaskFlowSingleMu.h new file mode 100644 index 00000000000..8a2dfaf6e7b --- /dev/null +++ b/PWG/muon/AliAnalysisTaskFlowSingleMu.h @@ -0,0 +1,80 @@ +#ifndef ALIANALYSISTASKFLOWSINGLEMU_H +#define ALIANALYSISTASKFLOWSINGLEMU_H + +/* $Id: AliAnalysisTaskFlowSingleMu.h 55545 2012-04-04 07:16:39Z pcrochet $ */ + +// +// AliAnalysisTaskFlowSingleMu +// Analysis task for flow of single muons in the spectrometer +// +// Author: Diego Stocco +// + +#include "AliVAnalysisMuon.h" +#include "TRandom3.h" + +class TObjArray; +class TString; +class TArrayD; +class TAxis; +class AliMuonTrackCuts; + +class AliAnalysisTaskFlowSingleMu : public AliVAnalysisMuon { + public: + AliAnalysisTaskFlowSingleMu(); + AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts); + virtual ~AliAnalysisTaskFlowSingleMu(); + + virtual void Terminate(Option_t *option); + + void MyUserCreateOutputObjects(); + void ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality); + + /* + enum { + kEPV0A, ///< EP form V0A + kEPTPC, ///< EP form TPC + kEPrandom ///< Random EP + }; + + void SetEPtype ( Int_t epType = kEPV0A ) { fEPtype = epType; } + */ + + private: + + AliAnalysisTaskFlowSingleMu(const AliAnalysisTaskFlowSingleMu&); + AliAnalysisTaskFlowSingleMu& operator=(const AliAnalysisTaskFlowSingleMu&); + + TArrayD GetCentralityRange(TString sRange); + + /* + enum { + kTrackContainer, ///< CF container for tracks + kHistoEP, ///< Event plane distribution + kNobjectTypes ///< Number of objects + }; + */ + + enum { + kStepReconstructed, ///< Reconstructed tracks + kStepGeneratedMC, ///< Generated tracks (MC) + kNsteps ///< Number of steps + }; + + enum { + kHvarPt, ///< Pt at vertex + kHvarEta, ///< Pseudo-Rapidity + kHvarPhi, ///< Phi + kHvarDeltaPhi, ///< Phi_mu - Psi_plane + kHvarCharge, ///< Particle charge + kHvarMotherType, ///< Mother type (MC only) + kNvars ///< THnSparse dimensions + }; + + TObjArray* fEPKeys; ///< EP keys + TRandom3* fRandom; //!< Random number generator + + ClassDef(AliAnalysisTaskFlowSingleMu, 1); // Single muon analysis +}; + +#endif -- 2.43.0