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