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