Fixes for std:: need with the trunk of Root (Jochen, Yves)
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskFlowSingleMu.cxx
CommitLineData
7e70f6e0 1/**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliAnalysisTaskFlowSingleMu.cxx 55545 2012-04-04 07:16:39Z pcrochet $ */
17
18//-----------------------------------------------------------------------------
19/// \class AliAnalysisTaskFlowSingleMu
20/// Analysis task for flow of single muons in the spectrometer.
21/// The output is a list of histograms and CF containers.
22/// The macro class can run on AODs or ESDs.
23/// If Monte Carlo information is present, some basics checks are performed.
24///
25/// \author Diego Stocco
26//-----------------------------------------------------------------------------
27
28#define AliAnalysisTaskFlowSingleMu_cxx
29
30#include "AliAnalysisTaskFlowSingleMu.h"
31
32// ROOT includes
33#include "TROOT.h"
34#include "TH1.h"
35#include "TH2.h"
36#include "TGraphErrors.h"
37#include "TAxis.h"
38#include "TCanvas.h"
39#include "TLegend.h"
40#include "TMath.h"
41#include "TObjString.h"
42#include "TObjArray.h"
43#include "TF1.h"
44#include "TStyle.h"
45//#include "TMCProcess.h"
46#include "TFitResultPtr.h"
47#include "TFile.h"
48#include "TArrayI.h"
49#include "TArrayD.h"
50#include "TSystem.h"
51#include "TGraphErrors.h"
52#include "TLatex.h"
53
54// STEER includes
55#include "AliAODEvent.h"
56#include "AliAODTrack.h"
57#include "AliAODMCParticle.h"
58#include "AliMCEvent.h"
59#include "AliMCParticle.h"
60#include "AliESDEvent.h"
61#include "AliESDMuonTrack.h"
62#include "AliVHeader.h"
63#include "AliAODMCHeader.h"
64#include "AliStack.h"
65#include "AliEventplane.h"
66
67// ANALYSIS includes
68#include "AliAnalysisManager.h"
69
70// CORRFW includes
71#include "AliCFContainer.h"
72#include "AliCFGridSparse.h"
73#include "AliCFEffGrid.h"
74
75// PWG includes
76#include "AliVAnalysisMuon.h"
77#include "AliMergeableCollection.h"
78#include "AliCounterCollection.h"
79#include "AliMuonTrackCuts.h"
80
81
82/// \cond CLASSIMP
83ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
84/// \endcond
85
e61afb80 86using std::ifstream;
7e70f6e0 87
88//________________________________________________________________________
89AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
90 AliVAnalysisMuon(),
91 fEPKeys(0x0),
92 fRandom(0x0)
93{
94 /// Default ctor.
95}
96
97//________________________________________________________________________
98AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
99 AliVAnalysisMuon(name, cuts),
100 fEPKeys(0x0),
101 fRandom(0x0)
102{
103 //
104 /// Constructor.
105 //
106 TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
107 fEPKeys = epMethods.Tokenize(" ");
108}
109
110
111//________________________________________________________________________
112AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
113{
114 //
115 /// Destructor
116 //
117 delete fEPKeys;
118 delete fRandom;
119}
120
121//___________________________________________________________________________
122void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
123{
124 //
125 /// Create outputs
126 //
127
128 if ( ! fRandom ) fRandom = new TRandom3();
129 // Set seed here, to avoid having the same seed for all jobs in grid
130 fRandom->SetSeed(0);
131
132 Int_t nPtBins = 160;
133 Double_t ptMin = 0., ptMax = 80.;
134 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
135
136 Int_t nEtaBins = 25;
137 Double_t etaMin = -4.5, etaMax = -2.;
138 TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
139
140 Int_t nPhiBins = 36;
141 Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
142 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
143
144 Int_t nDeltaPhiBins = 18;
145 Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
146 TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
147
148 Int_t nChargeBins = 2;
149 Double_t chargeMin = -2, chargeMax = 2.;
150 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
151
152 Int_t nMotherTypeBins = kNtrackSources;
153 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
154 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
155
156 Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
157 Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
158 Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
159 TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
160 TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
161
162 Int_t iobj=0;
163 TH1* histo = 0x0;
164 TString currTitle = "";
165 for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
166 TString currEP = fEPKeys->At(iep)->GetName();
167 AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
168
169 TString currName = "";
170 for ( Int_t idim = 0; idim<kNvars; idim++){
171 currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
172 currName.ReplaceAll("()","");
173
174 cfContainer->SetVarTitle(idim, currName.Data());
175 cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
176 }
177
178 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
179
180 TAxis* currAxis = 0x0;
181 for (Int_t istep=0; istep<kNsteps; istep++){
182 cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
183 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
184
185 currAxis = gridSparse->GetAxis(kHvarMotherType);
186 for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
187 currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
188 }
189 }
190
191 AddObjectToCollection(cfContainer, iobj++);
192
193 Bool_t isTPC = currEP.Contains("TPC");
194 Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
195 Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
196
197 histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
198 histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
199 AddObjectToCollection(histo, iobj++);
200
201 for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
202 TString auxEP = fEPKeys->At(jep)->GetName();
203 currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
204 histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
205 histo->GetXaxis()->SetTitle(currTitle.Data());
206 AddObjectToCollection(histo, iobj++);
207 }
208
209 currTitle = Form("cos(#psi_{1} - #psi_{2})");
210 histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
211 histo->GetXaxis()->SetTitle(currTitle.Data());
212
213 AddObjectToCollection(histo, iobj++);
214 } // loop on EPs
215
216 fMuonTrackCuts->Print("mask");
217}
218
219//________________________________________________________________________
220void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
221{
222 //
223 /// Fill output objects
224 //
225
226 //
227 // Base event cuts
228 //
229 if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
230 //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
231 //if ( isPileupFromSPD ) return;
232 if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) return;
233
234 TArrayD psiPlane(fEPKeys->GetEntriesFast());
235 psiPlane.Reset(-999.);
236 const Double_t kLowestVal = -10.;
237 Int_t harmonic = 2;
238
239 //
240 // Fill EP info
241 //
242 AliEventplane* eventPlane = InputEvent()->GetEventplane();
243 psiPlane[0] = eventPlane->GetEventplane("V0A",InputEvent(),harmonic);
244 psiPlane[1] = eventPlane->GetEventplane("V0C",InputEvent(),harmonic);
245 psiPlane[2] = eventPlane->GetEventplane("Q");
246 if ( psiPlane[2] < 0.) psiPlane[2] = -999.;
247 Double_t qx = 0., qy = 0.;
248 psiPlane[3] = eventPlane->CalculateVZEROEventPlane(InputEvent(),0,harmonic,qx,qy); // V0C ring 1
249 psiPlane[4] = eventPlane->CalculateVZEROEventPlane(InputEvent(),3,harmonic,qx,qy); // V0C ring 3
250 //psiPlane[3] = fRandom->Uniform(-TMath::Pi()/2., TMath::Pi()/2.);
251
252 Bool_t noValidEP = kTRUE;
253 TArrayD psiPlaneCalc(fEPKeys->GetEntriesFast());
254 for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
255 psiPlaneCalc[iep] = psiPlane[iep];
256 if ( psiPlane[iep] < kLowestVal ) continue;
257 if ( psiPlane[iep] < 0. ) psiPlaneCalc[iep] += TMath::Pi();
258 noValidEP = kFALSE;
259 }
260 if ( noValidEP ) return;
261
262 for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
263 if ( psiPlane[iep] < kLowestVal ) continue;
264 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
265 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
266 TString currEP = fEPKeys->At(iep)->GetName();
267 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hEventPlane%s",currEP.Data())))->Fill(psiPlane[iep]);
268 for ( Int_t jep=iep+1; jep<fEPKeys->GetEntriesFast(); jep++ ) {
269 if ( psiPlane[jep] < kLowestVal ) continue;
270 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub3%s%s",currEP.Data(), fEPKeys->At(jep)->GetName())))->Fill(TMath::Cos(2.*(psiPlaneCalc[iep]-psiPlaneCalc[jep])));
271 } // loop on auxialiary EP
272 if ( iep == 2 ) {
273 TVector2* qsub1 = eventPlane->GetQsub1();
274 TVector2* qsub2 = eventPlane->GetQsub2();
275 if ( qsub1 && qsub2 )
276 ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.));
277 }
278 } // loop on trigger
279 } // loop on EP
280
281 //
282 // Fill track info
283 //
284 Double_t containerInput[kNvars];
285 AliVParticle* track = 0x0;
286 for ( Int_t istep = 0; istep<2; ++istep ) {
287 Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
288 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
289 track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
290
291 Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
292 if ( ! isSelected ) continue;
293
294 Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
295
296
297 containerInput[kHvarPt] = track->Pt();
298 containerInput[kHvarEta] = track->Eta();
299 containerInput[kHvarPhi] = track->Phi();
300 containerInput[kHvarCharge] = track->Charge()/3.;
301 containerInput[kHvarMotherType] = (Double_t)trackSrc;
302
303 for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
304 if ( psiPlane[iep] < kLowestVal ) continue;
305 TString currEP = fEPKeys->At(iep)->GetName();
306 Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
307
308 // The difference between phi - psiPlane must be between (0, pi)
309 // Since the reaction plain is symmetric in pi,
310 // we can do in such a way that:
311 // -pi < delta_phi + N*pi < pi
312 // We choose N accrodingly
313 Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
314 deltaPhi += nPi * TMath::Pi();
315 containerInput[kHvarDeltaPhi] = deltaPhi;
316
317 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
318 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
319 if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
320 ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
321 } // loop on selected trigger classes
322 } // loop on event plane
323 } // loop on tracks
324 } // loop on container steps
325}
326
327
328//________________________________________________________________________
329TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
330{
331 //
332 /// Get the range from string
333 //
334 TArrayD centralityRange(2);
335 centralityRange.Reset(-999.);
336 TObjArray* array = sRange.Tokenize("_");
337 if ( array->GetEntries() == 2 ) {
338 for ( Int_t ientry=0; ientry<2; ientry++ ) {
339 centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
340 }
341 }
342 delete array;
343 return centralityRange;
344}
345
346//________________________________________________________________________
347void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
348 //
349 /// Draw some histograms at the end.
350 //
351
352 AliVAnalysisMuon::Terminate("");
353
354 if ( ! fMergeableCollection ) return;
355
356 TString physSel = fTerminateOptions->At(0)->GetName();
357 TString trigClassName = fTerminateOptions->At(1)->GetName();
358 TString centralityRange = fTerminateOptions->At(2)->GetName();
359 TString furtherOpt = fTerminateOptions->At(3)->GetName();
360
361 TObjArray* centralityClasses = centralityRange.Tokenize(",");
362 TArrayD centralityArray(centralityClasses->GetEntries()+1);
363 Int_t ib = 0;
364 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
365 TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
366 if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
367 centralityArray.AddAt(range[1],ib++);
368 }
369
370 TString selectEP = "";
371 TObjArray resoEParray;
372 resoEParray.SetOwner();
373
374 TString outFilename = "";
375 TObjArray* optArr = furtherOpt.Tokenize(" ");
376 TString currName = "";
377 TString resoFilename = "";
378 for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
379 currName = optArr->At(iopt)->GetName();
380 if ( currName.Contains(".root") ) outFilename = currName;
381 else if ( currName.Contains(".txt") ) resoFilename = currName;
382 else {
383 for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
384 if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
385 }
386 }
387 }
388 delete optArr;
389
390 if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
391
392 furtherOpt.ToUpper();
393
394 Bool_t v2only = furtherOpt.Contains("V2ONLY");
395
396 TAxis* customBins = 0x0;
397
398 //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
399 Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
400 //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
401 //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
402 //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
403 Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
404 customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
405
406 TCanvas *can = 0x0, *showCan = 0x0;
407 Int_t xshift = 100;
408 Int_t yshift = 20;
409 Int_t igroup1 = 0;
410 Int_t igroup2 = 0;
411
412 //
413 /// Read plane resolution file (if present)
414 //
415 TObjArray resoCentrality;
416 resoCentrality.SetOwner();
417 TArrayD resolution[3];
418 for ( Int_t ires=0; ires<3; ires++ ) {
419 resolution[ires].Set(50);
420 }
421
422 Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
423 Bool_t internalReso = kTRUE;
424 if ( externalReso ) {
425 // If plane resolution file exists, fill resolution
426 TString currLine = "";
427 ifstream inFile(resoFilename.Data());
428 if (inFile.is_open()) {
429 ib = 0;
430 while (! inFile.eof() ) {
431 currLine.ReadLine(inFile,kFALSE);
432 if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
433 TObjArray* array = currLine.Tokenize(" ");
434 resoCentrality.Add(new TObjString(array->At(0)->GetName()));
435 for ( Int_t ires=0; ires<3; ires++ ) {
436 resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
437 }
438 delete array;
439 ib++;
440 } // ! EOF
441 } // file is open
442 for ( Int_t ires=0; ires<3; ires++ ) {
443 resolution[ires].Set(ib);
444 }
445 } // file exists
446 else {
447 // If no plane resolution file exist,
448 // use the centralities in terminateOption
449 // and set the resolution to 1 and errors to 0.
450
451 for ( Int_t ires=0; ires<3; ires++ ) {
452 resolution[ires].Set(centralityClasses->GetEntries());
453 Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
454 resolution[ires].Reset(resetVal);
455 }
456 TArrayD currCos(3);
457 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
458 resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
459 Int_t icos = 0;
460 Double_t sumRelErrSquare = 0.;
461 for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
462 for ( Int_t jsub=0; jsub<resoEParray.GetEntries(); jsub++ ) {
463 TH1* histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName())));
464 if ( ! histo || histo->Integral() == 0 ) continue;
465 currCos[icos] = histo->GetMean();
466 if ( currCos[icos] <= 0. ) continue;
467 Double_t relErr = histo->GetMeanError() / currCos[icos];
468 sumRelErrSquare += relErr * relErr;
469 icos++;
470 } // loop on sub events
471 } // loop on sub events
472 if ( icos != 3 ) {
473 printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
474 internalReso = kFALSE;
475 continue;
476 }
477 resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
478 resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
479 printf("Resolution %s %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
480 }
481 }
482
483 Bool_t hasResolution = ( externalReso || internalReso );
484
485
486 TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
487 ib = 0;
488 for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
489 TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
490 if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
491 resoCentrArray.AddAt(range[1], ib++);
492 }
493
494 if ( hasResolution ) {
495 currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
496 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
497 TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries());
498 for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
499
500 resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
501 resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
502 }
503 resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
504 resoGraph->GetXaxis()->SetTitle("Centrality");
505 resoGraph->GetYaxis()->SetTitle("Resolution");
506 resoGraph->Draw("apz");
507 }
508
509 const Int_t kNresoCentr = resoCentrality.GetEntries();
510
511 //
512 // Create gridSparse array
513 //
514 AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
515 TString stepName[2];
516 for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
517 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
518 gridSparseArray[istep*kNresoCentr+icent] = 0x0;
519 }
520
521 AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
522 if ( ! cfContainer ) continue;
523 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
524 stepName[istep] = cfContainer->GetStepTitle(istep);
525 gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
526 } // loop on step
527 } // loop on centrality
528
529 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
530
531 Bool_t isMC = furtherOpt.Contains("MC");
532 Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
533 Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
534 if ( ! isMC ) srcColors[kUnidentified] = 1;
535
536 TList outList;
537
538 TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
539 Int_t drawOrder[3] = {0, 2, 1};
540 TString drawOpt[3] = {"apz","pz","a2"};
541 Int_t nTypes = ( hasResolution ) ? 3 : 1;
542
543 TString histoName = "", histoPattern = "";
544 ///////////
545 // Flow //
546 ///////////
547 gStyle->SetOptFit(1111);
548 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)))";
549 //
550 TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
551 if ( v2only ) func->SetParNames("scale", "v2");
552 else func->SetParNames("scale","v1", "v2", "v3");
553 Int_t v2par = ( v2only ) ? 1 : 2;
554
555 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
556 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
557 TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
558 TGraphErrors* flowVsCentrality[3];
559 for ( Int_t itype=0; itype<nTypes; itype++ ) {
560 histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
561 flowVsCentrality[itype] = new TGraphErrors();
562 flowVsCentrality[itype]->SetName(histoName.Data());
563 flowVsCentrality[itype]->SetTitle(histoName.Data());
564 }
565 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
566 TH2* hDeltaPhi2D = 0x0;
567 TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
568 TArrayD weightedReso(3);
569 weightedReso.Reset(0.);
570 Double_t nMuons = 0.;
571 TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
572 for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
573 if ( resoCentrArray[rcent] < centralityArray[icent] ||
574 resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
575 AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
576 if ( ! gridSparse || gridSparse->GetEntries() == 0. ) continue;
577
578 SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
579 SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
580 TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
581
582 Double_t currYield = auxHisto->Integral();
583 if ( currYield == 0. ) continue;
584 nMuons += currYield;
585 debugString += Form("\nAdding %s yield %g reso", resoCentrality.At(rcent)->GetName(), currYield);
586 for ( Int_t ires=0; ires<3; ires++ ) {
587 weightedReso[ires] += resolution[ires][rcent] * currYield;
588 debugString += Form(" %g", resolution[ires][rcent]);
589 }
590 histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data());
591 if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast<TH2*>(auxHisto->Clone(histoName.Data()));
592 else hDeltaPhi2D->Add(auxHisto);
593 delete auxHisto;
594 }
595
596 if ( ! hDeltaPhi2D ) continue;
597
598 debugString += Form("\nWeighted reso ");
599 for ( Int_t ires=0; ires<3; ires++ ) {
600 weightedReso[ires] = weightedReso[ires]/nMuons;
601 debugString += Form(" %g", weightedReso[ires]);
602 }
603 AliDebug(1,debugString.Data());
604
605 TAxis* ptAxis = hDeltaPhi2D->GetYaxis();
606
607 Int_t ipad = 0;
608 TGraphErrors* flowVsPt[3];
609 for ( Int_t itype=0; itype<nTypes; itype++ ) {
610 histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
611 flowVsPt[itype] = new TGraphErrors();
612 flowVsPt[itype]->SetName(histoName.Data());
613 flowVsPt[itype]->SetTitle(histoName.Data());
614 }
615
616 if ( ! customBins ) customBins = static_cast<TAxis*>(ptAxis->Clone());
617
618
619 for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
620 Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
621 Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
622 Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
623 Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
624 Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
625 Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
626 TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin);
627 projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
628 if ( projHisto->Integral() < 50 ) break;
629 if ( ipad%4 == 0 ) {
630 currName = projHisto->GetName();
631 currName.Append("_can");
632 showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
633 showCan->Divide(2,2);
634 ipad = 0;
635 }
636 showCan->cd(++ipad);
637 if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
638 else {
639 func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
640 func->FixParameter(1,0.);
641 }
642 projHisto->Fit(func,"R");
643 for ( Int_t itype=0; itype<nTypes; itype++ ) {
644 TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
645 Double_t resoVal = ( itype == 0 ) ? 1. : weightedReso[0];
646 Double_t resoErr = ( itype == 0 ) ? 0. : weightedReso[itype];
647 Double_t rawVal = func->GetParameter(v2par);
648 Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
649 Double_t yVal = rawVal/resoVal;
650 Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
651 Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
652 Int_t ipoint = currGraph->GetN();
653 currGraph->SetPoint(ipoint,xVal,yVal);
654 if ( itype > 0 && ipt != 0 ) {
655 // For the plot vs. pt the resolution error is fully correlated.
656 // Hence we do not use add it bin by bin, but rather write it
657 // on the title
658 resoErr = 0.;
659 TString graphTitle = currGraph->GetTitle();
660 if ( ! graphTitle.Contains("|") ) {
661 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]);
662 currGraph->SetTitle(graphTitle.Data());
663 }
664 }
665 Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
666 Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
667 Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
668 currGraph->SetPointError(ipoint,xErr,yErr);
669 } // loop on type
670 } // loop on pt bins
671 for ( Int_t itype=0; itype<nTypes; itype++ ) {
672 Int_t currType = drawOrder[itype];
673 if ( itype < 2 ) {
674 currName = flowVsPt[currType]->GetName();
675 currName.Append("Can");
676 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
677 igroup2++;
678 can->cd();
679 }
680 flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
681 flowVsPt[currType]->GetYaxis()->SetTitle("v2");
682 flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
683 flowVsPt[currType]->SetFillStyle(0);
684 flowVsPt[currType]->Draw(drawOpt[currType].Data());
685 outList.Add(flowVsPt[currType]);
686 } // loop on types
687 } // loop on centrality
688 for ( Int_t itype=0; itype<nTypes; itype++ ) {
689 Int_t currType = drawOrder[itype];
690 if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
691 if ( itype < 2 ) {
692 currName = flowVsCentrality[currType]->GetName();
693 currName.Append("Can");
694 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
695 can->cd();
696 }
697 flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
698 flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
699 flowVsCentrality[currType]->SetFillStyle(0);
700 flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
701 flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
702 outList.Add(flowVsCentrality[currType]);
703 igroup2++;
704 } // loop on type
705 igroup1++;
706 igroup2 = 0;
707 } // loop on track sources
708 } // loop on steps
709
710 delete customBins;
711
712
713 ///////////////////////
714 // Event plane check //
715 ///////////////////////
716 igroup1++;
717 igroup2 = 0;
718 Int_t icolor=0;
719 currName ="eventPlaneDistrib";
720 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
721 can->SetGridy();
722 igroup2++;
723
724 Bool_t addOffsetToCentrality = kFALSE;
725
726 TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
727 leg->SetBorderSize(1);
728 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
729 TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
730 if ( ! histo ) continue;
731 histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
732 if ( histo->Integral() < 50. ) continue;
733 if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
734 histo->Fit("pol0","R0");
735 TF1* pol0 = histo->GetFunction("pol0");
736 histo->Scale(1./pol0->GetParameter(0));
737 Int_t offset = ( addOffsetToCentrality ) ? icolor : 0;
738 TGraphErrors* graph = new TGraphErrors();
739 graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
740 for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
741 Int_t ibin = ipoint+1;
742 graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
743 graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
744 }
745 graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
746 graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
747 TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
748 if ( addOffsetToCentrality ) {
749 graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries());
750 graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
751 legTitle += Form(" ( + %i)", offset);
752 }
753 Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
754 graph->SetLineColor(currColor);
755 graph->SetMarkerColor(currColor);
756 graph->SetMarkerStyle(20+icolor);
757 TString currDrawOpt = "pz";
758 if ( icolor == 0 ) currDrawOpt += "a";
759 graph->Draw(currDrawOpt.Data());
760 legTitle.ReplaceAll("_"," - ");
761 leg->AddEntry(graph,legTitle.Data(),"lp");
762 if ( addOffsetToCentrality ) {
763 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()));
764 latex->SetTextColor(currColor);
765 latex->SetTextSize(0.025);
766 latex->Draw("same");
767 }
768 icolor++;
769 }
770 leg->Draw("same");
771
772
773 //////////////////////
774 // Event statistics //
775 //////////////////////
776 printf("\nTotal analyzed events:\n");
777 TString evtSel = Form("trigger:%s", trigClassName.Data());
778 fEventCounters->PrintSum(evtSel.Data());
779 printf("Physics selected analyzed events:\n");
780 evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
781 fEventCounters->PrintSum(evtSel.Data());
782
783 TString countPhysSel = "any";
784 if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
785 else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
786 countPhysSel.Prepend("selected:");
787 printf("Analyzed events vs. centrality:\n");
788 evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
789 fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
790
791 if ( ! outFilename.IsNull() ) {
792 printf("\nWriting output file %s\n", outFilename.Data());
793 TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
794 outList.Write();
795 file->Close();
796 }
797}