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