AliAODEvent::GetHeader now return AliVHeader
[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"
ebba9ef0 79#include "AliMuonEventCuts.h"
7e70f6e0 80#include "AliMuonTrackCuts.h"
ebba9ef0 81#include "AliAnalysisMuonUtility.h"
7e70f6e0 82
83
84/// \cond CLASSIMP
85ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
86/// \endcond
87
e61afb80 88using std::ifstream;
7e70f6e0 89
90//________________________________________________________________________
91AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
92 AliVAnalysisMuon(),
93 fEPKeys(0x0),
94 fRandom(0x0)
95{
96 /// Default ctor.
97}
98
99//________________________________________________________________________
100AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
101 AliVAnalysisMuon(name, cuts),
102 fEPKeys(0x0),
103 fRandom(0x0)
104{
105 //
106 /// Constructor.
107 //
108 TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
109 fEPKeys = epMethods.Tokenize(" ");
110}
111
112
113//________________________________________________________________________
114AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
115{
116 //
117 /// Destructor
118 //
119 delete fEPKeys;
120 delete fRandom;
121}
122
123//___________________________________________________________________________
124void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
125{
126 //
127 /// Create outputs
128 //
129
130 if ( ! fRandom ) fRandom = new TRandom3();
131 // Set seed here, to avoid having the same seed for all jobs in grid
132 fRandom->SetSeed(0);
133
134 Int_t nPtBins = 160;
135 Double_t ptMin = 0., ptMax = 80.;
136 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
137
138 Int_t nEtaBins = 25;
139 Double_t etaMin = -4.5, etaMax = -2.;
140 TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
141
142 Int_t nPhiBins = 36;
143 Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
144 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
145
146 Int_t nDeltaPhiBins = 18;
147 Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
148 TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
149
150 Int_t nChargeBins = 2;
151 Double_t chargeMin = -2, chargeMax = 2.;
152 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
153
154 Int_t nMotherTypeBins = kNtrackSources;
155 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
156 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
157
158 Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
159 Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
160 Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
161 TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
162 TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
163
164 Int_t iobj=0;
165 TH1* histo = 0x0;
166 TString currTitle = "";
167 for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
168 TString currEP = fEPKeys->At(iep)->GetName();
169 AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
170
171 TString currName = "";
172 for ( Int_t idim = 0; idim<kNvars; idim++){
173 currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
174 currName.ReplaceAll("()","");
175
176 cfContainer->SetVarTitle(idim, currName.Data());
177 cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
178 }
179
180 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
181
182 TAxis* currAxis = 0x0;
183 for (Int_t istep=0; istep<kNsteps; istep++){
184 cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
185 AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
186
187 currAxis = gridSparse->GetAxis(kHvarMotherType);
188 for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
189 currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
190 }
191 }
192
193 AddObjectToCollection(cfContainer, iobj++);
194
195 Bool_t isTPC = currEP.Contains("TPC");
196 Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
197 Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
198
199 histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
200 histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
201 AddObjectToCollection(histo, iobj++);
202
203 for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
204 TString auxEP = fEPKeys->At(jep)->GetName();
205 currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
206 histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
207 histo->GetXaxis()->SetTitle(currTitle.Data());
208 AddObjectToCollection(histo, iobj++);
209 }
210
211 currTitle = Form("cos(#psi_{1} - #psi_{2})");
212 histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
213 histo->GetXaxis()->SetTitle(currTitle.Data());
214
215 AddObjectToCollection(histo, iobj++);
216 } // loop on EPs
217
218 fMuonTrackCuts->Print("mask");
219}
220
221//________________________________________________________________________
222void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
223{
224 //
225 /// Fill output objects
226 //
227
228 //
229 // Base event cuts
230 //
7e70f6e0 231 //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
232 //if ( isPileupFromSPD ) return;
7e70f6e0 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 }
ebba9ef0 278 //else ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(2.*psiPlane[iep])); // REMEMBER TO CUT
7e70f6e0 279 } // loop on trigger
280 } // loop on EP
281
282 //
283 // Fill track info
284 //
285 Double_t containerInput[kNvars];
286 AliVParticle* track = 0x0;
ac2e110b 287 Int_t nSteps = MCEvent() ? 2 : 1;
288 for ( Int_t istep = 0; istep<nSteps; ++istep ) {
289 Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks();
7e70f6e0 290 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
ac2e110b 291 track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
7e70f6e0 292
293 Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
294 if ( ! isSelected ) continue;
295
33898693 296 Int_t trackSrc = GetParticleType(track);
7e70f6e0 297
298
299 containerInput[kHvarPt] = track->Pt();
300 containerInput[kHvarEta] = track->Eta();
301 containerInput[kHvarPhi] = track->Phi();
302 containerInput[kHvarCharge] = track->Charge()/3.;
303 containerInput[kHvarMotherType] = (Double_t)trackSrc;
304
305 for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
306 if ( psiPlane[iep] < kLowestVal ) continue;
307 TString currEP = fEPKeys->At(iep)->GetName();
308 Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
309
310 // The difference between phi - psiPlane must be between (0, pi)
311 // Since the reaction plain is symmetric in pi,
312 // we can do in such a way that:
313 // -pi < delta_phi + N*pi < pi
314 // We choose N accrodingly
315 Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
316 deltaPhi += nPi * TMath::Pi();
317 containerInput[kHvarDeltaPhi] = deltaPhi;
318
319 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
320 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
ebba9ef0 321 if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
7e70f6e0 322 ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
323 } // loop on selected trigger classes
324 } // loop on event plane
325 } // loop on tracks
326 } // loop on container steps
327}
328
329
330//________________________________________________________________________
331TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
332{
333 //
334 /// Get the range from string
335 //
336 TArrayD centralityRange(2);
337 centralityRange.Reset(-999.);
338 TObjArray* array = sRange.Tokenize("_");
339 if ( array->GetEntries() == 2 ) {
340 for ( Int_t ientry=0; ientry<2; ientry++ ) {
341 centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
342 }
343 }
344 delete array;
345 return centralityRange;
346}
347
348//________________________________________________________________________
349void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
350 //
351 /// Draw some histograms at the end.
352 //
353
354 AliVAnalysisMuon::Terminate("");
355
356 if ( ! fMergeableCollection ) return;
357
358 TString physSel = fTerminateOptions->At(0)->GetName();
359 TString trigClassName = fTerminateOptions->At(1)->GetName();
360 TString centralityRange = fTerminateOptions->At(2)->GetName();
361 TString furtherOpt = fTerminateOptions->At(3)->GetName();
362
363 TObjArray* centralityClasses = centralityRange.Tokenize(",");
364 TArrayD centralityArray(centralityClasses->GetEntries()+1);
365 Int_t ib = 0;
366 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
367 TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
368 if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
369 centralityArray.AddAt(range[1],ib++);
370 }
371
372 TString selectEP = "";
373 TObjArray resoEParray;
374 resoEParray.SetOwner();
375
376 TString outFilename = "";
377 TObjArray* optArr = furtherOpt.Tokenize(" ");
378 TString currName = "";
379 TString resoFilename = "";
380 for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
381 currName = optArr->At(iopt)->GetName();
382 if ( currName.Contains(".root") ) outFilename = currName;
383 else if ( currName.Contains(".txt") ) resoFilename = currName;
384 else {
385 for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
386 if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
387 }
388 }
389 }
390 delete optArr;
391
392 if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
393
394 furtherOpt.ToUpper();
395
396 Bool_t v2only = furtherOpt.Contains("V2ONLY");
397
398 TAxis* customBins = 0x0;
399
400 //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
ebba9ef0 401 //Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
402 //Double_t myBins[] = {6., 8., 10.};
403 Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
7e70f6e0 404 //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
405 //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
406 Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
407 customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
408
409 TCanvas *can = 0x0, *showCan = 0x0;
410 Int_t xshift = 100;
411 Int_t yshift = 20;
412 Int_t igroup1 = 0;
413 Int_t igroup2 = 0;
ebba9ef0 414
415 TList outList;
7e70f6e0 416
417 //
418 /// Read plane resolution file (if present)
419 //
420 TObjArray resoCentrality;
421 resoCentrality.SetOwner();
422 TArrayD resolution[3];
423 for ( Int_t ires=0; ires<3; ires++ ) {
424 resolution[ires].Set(50);
425 }
426
427 Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
428 Bool_t internalReso = kTRUE;
429 if ( externalReso ) {
430 // If plane resolution file exists, fill resolution
431 TString currLine = "";
2ad98fc5 432 std::ifstream inFile(resoFilename.Data());
7e70f6e0 433 if (inFile.is_open()) {
434 ib = 0;
435 while (! inFile.eof() ) {
436 currLine.ReadLine(inFile,kFALSE);
437 if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
438 TObjArray* array = currLine.Tokenize(" ");
439 resoCentrality.Add(new TObjString(array->At(0)->GetName()));
440 for ( Int_t ires=0; ires<3; ires++ ) {
441 resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
442 }
443 delete array;
444 ib++;
445 } // ! EOF
446 } // file is open
447 for ( Int_t ires=0; ires<3; ires++ ) {
448 resolution[ires].Set(ib);
449 }
450 } // file exists
451 else {
452 // If no plane resolution file exist,
453 // use the centralities in terminateOption
454 // and set the resolution to 1 and errors to 0.
455
456 for ( Int_t ires=0; ires<3; ires++ ) {
457 resolution[ires].Set(centralityClasses->GetEntries());
458 Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
459 resolution[ires].Reset(resetVal);
460 }
461 TArrayD currCos(3);
462 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
463 resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
464 Int_t icos = 0;
465 Double_t sumRelErrSquare = 0.;
ebba9ef0 466 TString hResoTitle = "";
7e70f6e0 467 for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
ebba9ef0 468 TString resoDet1 = resoEParray.At(isub)->GetName();
469 for ( Int_t jsub=isub+1; jsub<resoEParray.GetEntries(); jsub++ ) {
470 TString resoDet2 = resoEParray.At(jsub)->GetName();
471 TH1* histo = 0x0;
472 // Try both combinations: we want the first two values to contain the selectedEP!
473 TString resoDet = "";
474 for ( Int_t icomb=0; icomb<2; icomb++ ) {
475 TString hName = "hResoSub3";
476 resoDet = ( icomb == 0 ) ? resoDet1 + resoDet2 : resoDet2 + resoDet1;
477 hName += resoDet;
478 histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),hName));
479 if ( histo ) break;
480 }
7e70f6e0 481 if ( ! histo || histo->Integral() == 0 ) continue;
482 currCos[icos] = histo->GetMean();
483 if ( currCos[icos] <= 0. ) continue;
484 Double_t relErr = histo->GetMeanError() / currCos[icos];
485 sumRelErrSquare += relErr * relErr;
486 icos++;
ebba9ef0 487 hResoTitle += Form("%s ", resoDet.Data());
7e70f6e0 488 } // loop on sub events
489 } // loop on sub events
490 if ( icos != 3 ) {
491 printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
492 internalReso = kFALSE;
493 continue;
494 }
495 resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
496 resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
ebba9ef0 497 printf("Resolution (%s) %s %g +- %g\n", hResoTitle.Data(), centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
7e70f6e0 498 }
499 }
500
501 Bool_t hasResolution = ( externalReso || internalReso );
502
503
504 TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
505 ib = 0;
506 for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
507 TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
508 if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
509 resoCentrArray.AddAt(range[1], ib++);
510 }
ebba9ef0 511
7e70f6e0 512 const Int_t kNresoCentr = resoCentrality.GetEntries();
513
514 //
515 // Create gridSparse array
516 //
517 AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
518 TString stepName[2];
519 for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
520 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
521 gridSparseArray[istep*kNresoCentr+icent] = 0x0;
522 }
523
524 AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
525 if ( ! cfContainer ) continue;
526 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
527 stepName[istep] = cfContainer->GetStepTitle(istep);
528 gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
ebba9ef0 529 if ( ! customBins ) customBins = gridSparseArray[istep*kNresoCentr+icent]->GetAxis(kHvarPt);
7e70f6e0 530 } // loop on step
531 } // loop on centrality
ebba9ef0 532
7e70f6e0 533
534 Bool_t isMC = furtherOpt.Contains("MC");
535 Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
536 Int_t lastSrc = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
ebba9ef0 537 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
7e70f6e0 538 if ( ! isMC ) srcColors[kUnidentified] = 1;
ebba9ef0 539
540
7e70f6e0 541 //
ebba9ef0 542 // Get histograms
543 //
544 const Int_t kNcentralities = centralityClasses->GetEntries();
545 const Int_t kNhistos = kNsteps*kNtrackSources*kNcentralities;
546 TObjArray ptVsDeltaPhiList(kNhistos);
547 TObjArray ptVsPhiList(kNhistos);
548 TArrayD wgtReso[3];
549 for ( Int_t ires=0; ires<3; ires++ ) {
550 wgtReso[ires].Set(kNhistos);
551 wgtReso[ires].Reset(0.);
552 }
7e70f6e0 553 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
554 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
555 TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
7e70f6e0 556 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
ebba9ef0 557 Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
558 TH1* hPtVsDeltaPhi = 0x0;
559 TH1* hPtVsPhi = 0x0;
7e70f6e0 560 TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
7e70f6e0 561 Double_t nMuons = 0.;
562 TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
563 for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
564 if ( resoCentrArray[rcent] < centralityArray[icent] ||
ebba9ef0 565 resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
7e70f6e0 566 AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
567 if ( ! gridSparse || gridSparse->GetEntries() == 0. ) continue;
ebba9ef0 568
7e70f6e0 569 SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
570 SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
571 TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
ebba9ef0 572
7e70f6e0 573 Double_t currYield = auxHisto->Integral();
574 if ( currYield == 0. ) continue;
575 nMuons += currYield;
576 debugString += Form("\nAdding %s yield %g reso", resoCentrality.At(rcent)->GetName(), currYield);
577 for ( Int_t ires=0; ires<3; ires++ ) {
ebba9ef0 578 wgtReso[ires][idx] += resolution[ires][rcent] * currYield;
7e70f6e0 579 debugString += Form(" %g", resolution[ires][rcent]);
580 }
ebba9ef0 581 if ( ! hPtVsDeltaPhi ) hPtVsDeltaPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsDeltaPhi_%s", baseNameCent.Data())));
582 else hPtVsDeltaPhi->Add(auxHisto);
7e70f6e0 583 delete auxHisto;
ebba9ef0 584
585 auxHisto = gridSparse->Project(kHvarPhi,kHvarPt);
586 if ( ! hPtVsPhi ) hPtVsPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsPhi_%s", baseNameCent.Data())));
587 else hPtVsPhi->Add(auxHisto);
588 delete auxHisto;
589
590 } // loop on resolution centralities
591
592 if ( nMuons == 0. ) continue;
7e70f6e0 593 debugString += Form("\nWeighted reso ");
594 for ( Int_t ires=0; ires<3; ires++ ) {
ebba9ef0 595 wgtReso[ires][idx] = wgtReso[ires][idx]/nMuons;
596 debugString += Form(" %g", wgtReso[ires][idx]);
7e70f6e0 597 }
598 AliDebug(1,debugString.Data());
ebba9ef0 599
600 ptVsDeltaPhiList.AddAt(hPtVsDeltaPhi,idx);
601 ptVsPhiList.AddAt(hPtVsPhi,idx);
602 } // loop on centralities
603 } // loop on sources
604 } // loop on steps
605
606
607
608 /////////////////////
609 // Plot resolution //
610 /////////////////////
611 if ( hasResolution ) {
612 currName = Form("reso_%s", selectEP.Data());
613 currName += externalReso ? Form("_external") : Form("_sub_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
614 TGraphErrors* resoGraph = new TGraphErrors(kNresoCentr);
615 resoGraph->SetName(currName.Data());
616 currName.Append("_can");
617 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
3674921d 618 can->cd();
ebba9ef0 619
620 for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
621 resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
622 resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
623 }
624 resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
625 resoGraph->GetXaxis()->SetTitle("Centrality");
626 resoGraph->GetYaxis()->SetTitle("Resolution");
627 resoGraph->Draw("apz");
628 outList.Add(resoGraph);
629
630 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
631 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
632 currName = Form("wgtReso_%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
633 resoGraph = new TGraphErrors();
634 resoGraph->SetName(currName.Data());
635 for ( Int_t icent=0; icent<kNcentralities; icent++ ) {
636 Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
637 if ( wgtReso[0][idx] == 0. ) continue;
638 resoGraph->SetPoint(icent, 0.5*(centralityArray[icent+1] + centralityArray[icent]), wgtReso[0][idx]);
639 resoGraph->SetPointError(icent, 0.5*(centralityArray[icent+1] - centralityArray[icent]), wgtReso[1][idx]);
640 } // loop on centralities
641 if ( resoGraph->GetN() == 0 ) continue;
642 printf("New reso %s\n", resoGraph->GetName());
643 resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
644 resoGraph->GetXaxis()->SetTitle("Centrality");
645 resoGraph->GetYaxis()->SetTitle("Resolution");
646 resoGraph->Draw("pz");
647 outList.Add(resoGraph);
648 } // loop on sources
649 } // loop on steps
650 }
651
652
653
654 TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
655 Int_t drawOrder[3] = {0, 2, 1};
656 TString drawOpt[3] = {"apz","pz","a2"};
657 Int_t nTypes = ( hasResolution ) ? 3 : 1;
7e70f6e0 658
ebba9ef0 659 TString histoName = "", histoPattern = "";
660 ///////////
661 // Flow //
662 ///////////
663
664 gStyle->SetOptFit(1111);
665 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)))";
666 //
667 TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
668 if ( v2only ) func->SetParNames("scale", "v2");
669 else func->SetParNames("scale","v1", "v2", "v3");
670 Int_t v2par = ( v2only ) ? 1 : 2;
671
672 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
673 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
674 TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
675 TGraphErrors* flowVsCentrality[3];
676 for ( Int_t itype=0; itype<nTypes; itype++ ) {
677 histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
678 flowVsCentrality[itype] = new TGraphErrors();
679 flowVsCentrality[itype]->SetName(histoName.Data());
680 flowVsCentrality[itype]->SetTitle(histoName.Data());
681 }
682 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
683 Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
684 TH2* histo2D = static_cast<TH2*> (ptVsDeltaPhiList.At(idx));
685 if ( ! histo2D ) continue;
686 TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
7e70f6e0 687
ebba9ef0 688 TAxis* ptAxis = histo2D->GetYaxis();
689
7e70f6e0 690 Int_t ipad = 0;
691 TGraphErrors* flowVsPt[3];
692 for ( Int_t itype=0; itype<nTypes; itype++ ) {
693 histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
694 flowVsPt[itype] = new TGraphErrors();
695 flowVsPt[itype]->SetName(histoName.Data());
696 flowVsPt[itype]->SetTitle(histoName.Data());
697 }
ebba9ef0 698
7e70f6e0 699 for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
700 Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
701 Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
702 Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
703 Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
704 Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
705 Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
ebba9ef0 706 TH1* projHisto = histo2D->ProjectionX(Form("%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
7e70f6e0 707 projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
708 if ( projHisto->Integral() < 50 ) break;
709 if ( ipad%4 == 0 ) {
710 currName = projHisto->GetName();
711 currName.Append("_can");
712 showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
713 showCan->Divide(2,2);
714 ipad = 0;
715 }
716 showCan->cd(++ipad);
717 if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
718 else {
719 func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
720 func->FixParameter(1,0.);
721 }
722 projHisto->Fit(func,"R");
723 for ( Int_t itype=0; itype<nTypes; itype++ ) {
724 TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
ebba9ef0 725 Double_t resoVal = ( itype == 0 ) ? 1. : wgtReso[0][idx];
726 Double_t resoErr = ( itype == 0 ) ? 0. : wgtReso[itype][idx];
7e70f6e0 727 Double_t rawVal = func->GetParameter(v2par);
728 Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
729 Double_t yVal = rawVal/resoVal;
730 Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
731 Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
732 Int_t ipoint = currGraph->GetN();
733 currGraph->SetPoint(ipoint,xVal,yVal);
734 if ( itype > 0 && ipt != 0 ) {
735 // For the plot vs. pt the resolution error is fully correlated.
736 // Hence we do not use add it bin by bin, but rather write it
737 // on the title
738 resoErr = 0.;
739 TString graphTitle = currGraph->GetTitle();
740 if ( ! graphTitle.Contains("|") ) {
ebba9ef0 741 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", wgtReso[1][idx]/wgtReso[0][idx], wgtReso[2][idx]/wgtReso[0][idx]);
7e70f6e0 742 currGraph->SetTitle(graphTitle.Data());
743 }
744 }
745 Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
746 Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
747 Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
748 currGraph->SetPointError(ipoint,xErr,yErr);
749 } // loop on type
750 } // loop on pt bins
751 for ( Int_t itype=0; itype<nTypes; itype++ ) {
752 Int_t currType = drawOrder[itype];
753 if ( itype < 2 ) {
754 currName = flowVsPt[currType]->GetName();
755 currName.Append("Can");
756 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
757 igroup2++;
758 can->cd();
759 }
760 flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
761 flowVsPt[currType]->GetYaxis()->SetTitle("v2");
762 flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
763 flowVsPt[currType]->SetFillStyle(0);
764 flowVsPt[currType]->Draw(drawOpt[currType].Data());
765 outList.Add(flowVsPt[currType]);
766 } // loop on types
767 } // loop on centrality
768 for ( Int_t itype=0; itype<nTypes; itype++ ) {
769 Int_t currType = drawOrder[itype];
770 if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
771 if ( itype < 2 ) {
772 currName = flowVsCentrality[currType]->GetName();
773 currName.Append("Can");
774 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
775 can->cd();
776 }
777 flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
778 flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
779 flowVsCentrality[currType]->SetFillStyle(0);
780 flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
781 flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
782 outList.Add(flowVsCentrality[currType]);
783 igroup2++;
784 } // loop on type
785 igroup1++;
786 igroup2 = 0;
ebba9ef0 787 } // loop on track sources
7e70f6e0 788 } // loop on steps
789
7e70f6e0 790 ///////////////////////
791 // Event plane check //
792 ///////////////////////
793 igroup1++;
794 igroup2 = 0;
795 Int_t icolor=0;
796 currName ="eventPlaneDistrib";
797 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
798 can->SetGridy();
799 igroup2++;
800
ebba9ef0 801 const Int_t kNfits = 2;
802 TString fitTypeName[kNfits] = {"Cos","Sin"};
803
7e70f6e0 804 Bool_t addOffsetToCentrality = kFALSE;
ebba9ef0 805 Double_t offsetStep = 0.1;
7e70f6e0 806
807 TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
808 leg->SetBorderSize(1);
ebba9ef0 809
810 TGraphErrors* epCheck[kNfits];
811 for ( Int_t itype=0; itype<kNfits; itype++ ) {
812 histoName = Form("checkFourier_%s",fitTypeName[itype].Data());
813 epCheck[itype] = new TGraphErrors();
814 epCheck[itype]->SetName(histoName.Data());
815// epCheck[itype]->SetTitle(histoName.Data());
816 }
7e70f6e0 817 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
818 TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
819 if ( ! histo ) continue;
820 histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
821 if ( histo->Integral() < 50. ) continue;
822 if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
ebba9ef0 823
824 for ( Int_t itype=0; itype<kNfits; itype++ ) {
825 TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
826 TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), histo->GetXaxis()->GetBinLowEdge(1), histo->GetXaxis()->GetBinUpEdge(histo->GetXaxis()->GetNbins()));
827 checkFunc->SetParameters(histo->Integral(), 0.);
828 histo->Fit(checkFunc,"R0");
829 TGraphErrors* currGraph = epCheck[itype];
830 Double_t yVal = checkFunc->GetParameter(1);
831 Double_t yErr = checkFunc->GetParError(1);
832 Double_t xVal = 0.5*(centralityArray[icent+1] + centralityArray[icent]);
833 Double_t xErr = 0.5*(centralityArray[icent+1] - centralityArray[icent]);
834 Int_t ipoint = currGraph->GetN();
835 currGraph->SetPoint(ipoint,xVal,yVal);
836 currGraph->SetPointError(ipoint,xErr,yErr);
837 } // loop on type
7e70f6e0 838 histo->Fit("pol0","R0");
839 TF1* pol0 = histo->GetFunction("pol0");
840 histo->Scale(1./pol0->GetParameter(0));
ebba9ef0 841 Double_t offset = ( addOffsetToCentrality ) ? offsetStep*(Double_t)icolor : 0.;
7e70f6e0 842 TGraphErrors* graph = new TGraphErrors();
843 graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
844 for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
845 Int_t ibin = ipoint+1;
846 graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
847 graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
848 }
849 graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
850 graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
851 TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
852 if ( addOffsetToCentrality ) {
ebba9ef0 853 graph->GetYaxis()->SetRangeUser(1.-offsetStep, 1.+1.5*offsetStep*(Double_t)centralityClasses->GetEntries());
7e70f6e0 854 graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
ebba9ef0 855 legTitle += Form(" ( + %g)", offset);
7e70f6e0 856 }
857 Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
858 graph->SetLineColor(currColor);
859 graph->SetMarkerColor(currColor);
860 graph->SetMarkerStyle(20+icolor);
861 TString currDrawOpt = "pz";
862 if ( icolor == 0 ) currDrawOpt += "a";
863 graph->Draw(currDrawOpt.Data());
864 legTitle.ReplaceAll("_"," - ");
865 leg->AddEntry(graph,legTitle.Data(),"lp");
866 if ( addOffsetToCentrality ) {
867 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()));
868 latex->SetTextColor(currColor);
869 latex->SetTextSize(0.025);
870 latex->Draw("same");
871 }
872 icolor++;
873 }
874 leg->Draw("same");
ebba9ef0 875
876 currName = "epCheckCan";
877 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
878 can->SetGridy();
879 igroup2++;
880 leg = new TLegend(0.7,0.7,0.95,0.95);
881 leg->SetBorderSize(1);
882 for ( Int_t itype=0; itype<kNfits; itype++ ) {
883 epCheck[itype]->SetLineColor(srcColors[itype]);
884 epCheck[itype]->SetMarkerColor(srcColors[itype]);
885 epCheck[itype]->SetMarkerStyle(20+itype);
886 epCheck[itype]->GetXaxis()->SetTitle("Centrality (%)");
887 TString currDrawOpt = "pz";
888 if ( can->GetListOfPrimitives()->GetEntries() == 0 ) currDrawOpt.Prepend("a");
889 epCheck[itype]->Draw(currDrawOpt.Data());
890 leg->AddEntry(epCheck[itype],Form("<%s(2 #psi_{%s})>", fitTypeName[itype].Data(), selectEP.Data()),"lp");
891 }
892 leg->Draw("same");
893
894 for ( Int_t istep=0; istep<kNsteps; istep++ ) {
895 for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
896 TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
897 for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
898 Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
899 TH2* histo2D = static_cast<TH2*> (ptVsPhiList.At(idx));
900 if ( ! histo2D ) continue;
901 TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
902
903 TAxis* ptAxis = histo2D->GetYaxis();
904
905 Int_t ipad = 0;
906 TGraphErrors* detEffect[kNfits];
907 for ( Int_t itype=0; itype<kNfits; itype++ ) {
908 histoName = Form("checkResoVsPt_%s_%s", baseNameCent.Data(), fitTypeName[itype].Data());
909 detEffect[itype] = new TGraphErrors();
910 detEffect[itype]->SetName(histoName.Data());
911 detEffect[itype]->SetTitle(histoName.Data());
912 }
913
914 for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
915 Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
916 Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
917 Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
918 Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
919 Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
920 Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
921 TH1* projHisto = histo2D->ProjectionX(Form("checkReso_%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
922 projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
923 if ( projHisto->Integral() < 50 ) break;
924 if ( ipad%4 == 0 ) {
925 currName = projHisto->GetName();
926 currName.Append("_can");
927 showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
928 showCan->Divide(2,2);
929 ipad = 0;
930 }
931 showCan->cd(++ipad);
932 for ( Int_t itype=0; itype<kNfits; itype++ ) {
933 TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
934 TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), projHisto->GetXaxis()->GetBinLowEdge(1), projHisto->GetXaxis()->GetBinUpEdge(projHisto->GetXaxis()->GetNbins()));
935 checkFunc->SetParameters(projHisto->Integral(), 0.);
936 projHisto->Fit(checkFunc,"R");
937 TGraphErrors* currGraph = detEffect[itype];
938 Double_t yVal = checkFunc->GetParameter(1);
939 Double_t yErr = checkFunc->GetParError(1);
940 Double_t xVal = 0.5*(ptMax + ptMin);
941 Double_t xErr = 0.5*(ptMax - ptMin);
942 Int_t ipoint = currGraph->GetN();
943 currGraph->SetPoint(ipoint,xVal,yVal);
944 currGraph->SetPointError(ipoint,xErr,yErr);
945 } // loop on type
946 } // loop on pt bins
947 for ( Int_t itype=0; itype<kNfits; itype++ ) {
948 currName = detEffect[itype]->GetName();
949 currName.Append("Can");
950 can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
951 igroup2++;
952 can->cd();
953 detEffect[itype]->GetXaxis()->SetTitle(ptAxis->GetTitle());
954 detEffect[itype]->GetYaxis()->SetTitle("v2");
955// detEffect[itype]->GetYaxis()->SetRangeUser(0., 0.25);
956 detEffect[itype]->SetFillStyle(0);
957 detEffect[itype]->Draw("ap");
958// outList.Add(detEffect[itype]);
959 } // loop on types
960 } // loop on centrality
961 } // loop on track sources
962 } // loop on steps
7e70f6e0 963
09d5920f 964 delete centralityClasses;
7e70f6e0 965
966 //////////////////////
967 // Event statistics //
968 //////////////////////
969 printf("\nTotal analyzed events:\n");
970 TString evtSel = Form("trigger:%s", trigClassName.Data());
971 fEventCounters->PrintSum(evtSel.Data());
972 printf("Physics selected analyzed events:\n");
973 evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
974 fEventCounters->PrintSum(evtSel.Data());
975
976 TString countPhysSel = "any";
977 if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
978 else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
979 countPhysSel.Prepend("selected:");
980 printf("Analyzed events vs. centrality:\n");
981 evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
982 fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
983
984 if ( ! outFilename.IsNull() ) {
985 printf("\nWriting output file %s\n", outFilename.Data());
986 TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
987 outList.Write();
988 file->Close();
989 }
ebba9ef0 990
991 delete customBins;
7e70f6e0 992}