]>
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; | |
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 | //________________________________________________________________________ | |
331 | TArrayD 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 | //________________________________________________________________________ | |
349 | void 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 | } |