]>
Commit | Line | Data |
---|---|---|
ce7adfe9 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2008, 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 | //***************************************************** | |
17 | // Class AliEPSelectionTask | |
18 | // Class to determine event plane | |
19 | // author: Alberica Toia, Johanna Gramling | |
20 | //***************************************************** | |
21 | ||
22 | #include "AliEPSelectionTask.h" | |
23 | ||
24 | #include <TTree.h> | |
25 | #include <TList.h> | |
26 | #include <TH1F.h> | |
27 | #include <TH2F.h> | |
28 | #include <TProfile.h> | |
29 | #include <TFile.h> | |
30 | #include <TObjString.h> | |
31 | #include <TString.h> | |
32 | #include <TCanvas.h> | |
33 | #include <TROOT.h> | |
34 | #include <TDirectory.h> | |
35 | #include <TSystem.h> | |
36 | #include <iostream> | |
37 | #include <TRandom2.h> | |
38 | #include <TArrayF.h> | |
39 | ||
40 | #include "AliAnalysisManager.h" | |
41 | #include "AliVEvent.h" | |
42 | #include "AliESD.h" | |
43 | #include "AliESDEvent.h" | |
44 | #include "AliMCEvent.h" | |
45 | #include "AliESDtrack.h" | |
46 | #include "AliESDtrackCuts.h" | |
47 | #include "AliESDHeader.h" | |
48 | #include "AliESDInputHandler.h" | |
49 | #include "AliAODHandler.h" | |
50 | #include "AliAODEvent.h" | |
51 | #include "AliHeader.h" | |
52 | #include "AliGenHijingEventHeader.h" | |
53 | #include "AliAnalysisTaskSE.h" | |
54 | #include "AliPhysicsSelectionTask.h" | |
55 | #include "AliPhysicsSelection.h" | |
56 | #include "AliBackgroundSelection.h" | |
57 | #include "AliESDUtils.h" | |
58 | ||
59 | #include "AliEventplane.h" | |
60 | ||
61 | ClassImp(AliEPSelectionTask) | |
62 | ||
63 | //________________________________________________________________________ | |
64 | AliEPSelectionTask::AliEPSelectionTask(): | |
65 | AliAnalysisTaskSE(), | |
66 | fDebug(0), | |
67 | fAnalysisInput("ESD"), | |
68 | fStatus("TPC"), | |
69 | fUseMCRP(kFALSE), | |
70 | fUsePhiWeight(kFALSE), | |
71 | fUsePtWeight(kFALSE), | |
72 | fSaveTrackContribution(kFALSE), | |
73 | fESDtrackCuts(0), | |
74 | ftracklist(0), | |
75 | fPhiDist(0), | |
76 | fQVector(0), | |
77 | fQContributionX(0), | |
78 | fQContributionY(0), | |
79 | fEventplaneQ(0), | |
80 | fQsub1(0), | |
81 | fQsub2(0), | |
82 | fQsubRes(0), | |
83 | fOutputList(0), | |
84 | fHOutEventplaneQ(0), | |
85 | fHOutPhi(0), | |
86 | fHOutPhiCorr(0), | |
87 | fHOutsub1sub2(0), | |
88 | fHOutNTEPRes(0), | |
89 | fHOutPTPsi(0), | |
90 | fHOutDiff(0), | |
91 | fHOutleadPTPsi(0) | |
92 | { | |
93 | // Default constructor | |
94 | AliInfo("Event Plane Selection enabled."); | |
95 | } | |
96 | ||
97 | //________________________________________________________________________ | |
98 | AliEPSelectionTask::AliEPSelectionTask(const char *name): | |
99 | AliAnalysisTaskSE(name), | |
100 | fDebug(0), | |
101 | fAnalysisInput("ESD"), | |
102 | fStatus("TPC"), | |
103 | fUseMCRP(kFALSE), | |
104 | fUsePhiWeight(kFALSE), | |
105 | fUsePtWeight(kFALSE), | |
106 | fSaveTrackContribution(kFALSE), | |
107 | fESDtrackCuts(0), | |
108 | ftracklist(0), | |
109 | fPhiDist(0), | |
110 | fQVector(0), | |
111 | fQContributionX(0), | |
112 | fQContributionY(0), | |
113 | fEventplaneQ(0), | |
114 | fQsub1(0), | |
115 | fQsub2(0), | |
116 | fQsubRes(0), | |
117 | fOutputList(0), | |
118 | fHOutEventplaneQ(0), | |
119 | fHOutPhi(0), | |
120 | fHOutPhiCorr(0), | |
121 | fHOutsub1sub2(0), | |
122 | fHOutNTEPRes(0), | |
123 | fHOutPTPsi(0), | |
124 | fHOutDiff(0), | |
125 | fHOutleadPTPsi(0) | |
126 | { | |
127 | // Default constructor | |
128 | AliInfo("Event Plane Selection enabled."); | |
129 | DefineOutput(1, TList::Class()); | |
130 | } | |
131 | ||
132 | //________________________________________________________________________ | |
133 | AliEPSelectionTask::~AliEPSelectionTask() | |
134 | { | |
135 | // Destructor | |
136 | if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){ | |
137 | delete fOutputList; | |
138 | fOutputList = 0; | |
139 | } | |
140 | if (fESDtrackCuts){ | |
141 | delete fESDtrackCuts; | |
142 | fESDtrackCuts = 0; | |
143 | } | |
144 | if (fQVector){ | |
145 | delete fQVector; | |
146 | fQVector = 0; | |
147 | } | |
148 | if (fQsub1){ | |
149 | delete fQsub1; | |
150 | fQsub1 = 0; | |
151 | } | |
152 | if (fQsub2){ | |
153 | delete fQsub2; | |
154 | fQsub2 = 0; | |
155 | } | |
156 | } | |
157 | ||
158 | //________________________________________________________________________ | |
159 | void AliEPSelectionTask::UserCreateOutputObjects() | |
160 | { | |
161 | // Create the output containers | |
162 | if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n"); | |
163 | AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo); | |
164 | ||
165 | fOutputList = new TList(); | |
166 | fOutputList->SetOwner(); | |
167 | fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi()); | |
168 | fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi()); | |
169 | fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi()); | |
170 | fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi()); | |
171 | fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi()); | |
172 | fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi()); | |
173 | fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi()); | |
174 | fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi()); | |
175 | ||
176 | fOutputList->Add(fHOutEventplaneQ); | |
177 | fOutputList->Add(fHOutPhi); | |
178 | fOutputList->Add(fHOutPhiCorr); | |
179 | fOutputList->Add(fHOutsub1sub2); | |
180 | fOutputList->Add(fHOutNTEPRes); | |
181 | fOutputList->Add(fHOutPTPsi); | |
182 | fOutputList->Add(fHOutleadPTPsi); | |
183 | fOutputList->Add(fHOutDiff); | |
184 | ||
185 | PostData(1, fOutputList); | |
186 | } | |
187 | ||
188 | //________________________________________________________________________ | |
189 | void AliEPSelectionTask::UserExec(Option_t */*option*/) | |
190 | { | |
191 | // Execute analysis for current event: | |
192 | if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n"); | |
193 | ||
194 | AliEventplane* esdEP = 0; | |
195 | TVector2 QQ1; | |
196 | TVector2 QQ2; | |
197 | Double_t fRP = 0.; // the monte carlo reaction plane angle | |
198 | ||
199 | if (fAnalysisInput.CompareTo("ESD")==0){ | |
200 | ||
201 | AliVEvent* event = InputEvent(); | |
202 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event); | |
203 | ||
204 | if (fUseMCRP) { | |
205 | AliMCEvent* mcEvent = MCEvent(); | |
206 | if (mcEvent && mcEvent->GenEventHeader()) { | |
207 | AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader()); | |
208 | fRP = headerH->ReactionPlaneAngle(); | |
209 | } | |
210 | } | |
211 | if (esd){ | |
212 | esdEP = esd->GetEventplane(); | |
213 | if (fSaveTrackContribution) { | |
214 | esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks()); | |
215 | esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks()); | |
216 | } | |
217 | ||
218 | if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE); | |
219 | if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE); | |
220 | const int NT = ftracklist->GetEntries(); | |
221 | ||
222 | if (NT>4){ | |
223 | fQVector = new TVector2(GetQ(esdEP,ftracklist)); | |
224 | fEventplaneQ = fQVector->Phi()/2; | |
225 | GetQsub(QQ1, QQ2, ftracklist); | |
226 | fQsub1 = new TVector2(QQ1); | |
227 | fQsub2 = new TVector2(QQ2); | |
228 | fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2); | |
ce7adfe9 | 229 | esdEP->SetQVector(fQVector); |
230 | esdEP->SetEventplaneQ(fEventplaneQ); | |
231 | esdEP->SetQsub(fQsub1,fQsub2); | |
232 | esdEP->SetQsubRes(fQsubRes); | |
233 | ||
234 | fHOutEventplaneQ->Fill(fEventplaneQ); | |
235 | fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2); | |
236 | fHOutNTEPRes->Fill(NT,fQsubRes); | |
237 | ||
238 | if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP); | |
239 | ||
240 | for (int iter = 0; iter<NT;iter++){ | |
241 | AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter)); | |
242 | float delta = track->Phi()-fEventplaneQ; | |
243 | while (delta < 0) delta += TMath::Pi(); | |
244 | while (delta > TMath::Pi()) delta -= TMath::Pi(); | |
245 | fHOutPTPsi->Fill(track->Pt(),delta); | |
246 | fHOutPhi->Fill(track->Phi()); | |
247 | fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track)); | |
248 | } | |
249 | ||
250 | AliESDtrack* trmax = esd->GetTrack(0); | |
251 | for (int iter = 1; iter<NT;iter++){ | |
252 | AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter)); | |
253 | if (track->Pt() > trmax->Pt()) trmax = track; | |
254 | } | |
255 | fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ); | |
256 | } | |
257 | delete ftracklist; | |
258 | } | |
259 | } | |
260 | ||
261 | else if (fAnalysisInput.CompareTo("AOD")==0){ | |
262 | //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent()); | |
263 | // to be implemented | |
264 | printf(" AOD analysis not yet implemented!!!\n\n"); | |
265 | return; | |
266 | } | |
267 | else { | |
268 | printf(" Analysis Input not known!\n\n "); | |
269 | return; | |
270 | } | |
271 | PostData(1, fOutputList); | |
272 | } | |
273 | ||
274 | //________________________________________________________________________ | |
275 | void AliEPSelectionTask::Terminate(Option_t */*option*/) | |
276 | { | |
277 | // Terminate analysis | |
278 | } | |
279 | ||
280 | //__________________________________________________________________________ | |
281 | TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist) | |
282 | { | |
283 | TVector2 mQ; | |
284 | float mQx=0, mQy=0; | |
285 | AliESDtrack* track; | |
286 | Double_t weight; | |
287 | ||
288 | int NT = tracklist->GetEntries(); | |
289 | ||
290 | for (int i=0; i<NT; i++){ | |
291 | weight = 1; | |
292 | track = dynamic_cast<AliESDtrack*> (tracklist->At(i)); | |
293 | weight = GetWeight(track); | |
294 | if (fSaveTrackContribution){ | |
295 | EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID()); | |
296 | EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID()); | |
297 | } | |
298 | mQx += (weight*cos(2*track->Phi())); | |
299 | mQy += (weight*sin(2*track->Phi())); | |
300 | } | |
301 | mQ.Set(mQx,mQy); | |
302 | return mQ; | |
303 | } | |
304 | ||
305 | //________________________________________________________________________ | |
306 | void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist) | |
307 | { | |
308 | TVector2 mQ[2]; | |
309 | float mQx1=0, mQy1=0, mQx2=0, mQy2=0; | |
310 | Double_t weight; | |
311 | ||
312 | AliESDtrack* track; | |
313 | TRandom2 RN = 0; | |
314 | ||
315 | int NT = tracklist->GetEntries(); | |
316 | int trackcounter1=0, trackcounter2=0; | |
317 | ||
318 | for (Int_t i = 0; i < NT; i++) { | |
319 | weight = 1; | |
320 | track = dynamic_cast<AliESDtrack*> (tracklist->At(i)); | |
321 | weight = GetWeight(track); | |
322 | ||
323 | // This loop splits the track set into 2 random subsets | |
324 | if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){ | |
325 | float random = RN.Rndm(); | |
326 | if(random < .5){ | |
327 | mQx1 += (weight*cos(2*track->Phi())); | |
328 | mQy1 += (weight*sin(2*track->Phi())); | |
329 | trackcounter1++; | |
330 | } | |
331 | else { | |
332 | mQx2 += (weight*cos(2*track->Phi())); | |
333 | mQy2 += (weight*sin(2*track->Phi())); | |
334 | trackcounter2++; | |
335 | } | |
336 | } | |
337 | else if( trackcounter1 >= int(NT/2.)){ | |
338 | mQx2 += (weight*cos(2*track->Phi())); | |
339 | mQy2 += (weight*sin(2*track->Phi())); | |
340 | trackcounter2++; | |
341 | } | |
342 | else { | |
343 | mQx1 += (weight*cos(2*track->Phi())); | |
344 | mQy1 += (weight*sin(2*track->Phi())); | |
345 | trackcounter1++; | |
346 | } | |
347 | } | |
348 | mQ[0].Set(mQx1,mQy1); | |
349 | mQ[1].Set(mQx2,mQy2); | |
350 | Q1 = mQ[0]; | |
351 | Q2 = mQ[1]; | |
352 | } | |
353 | ||
354 | //________________________________________________________________________ | |
355 | void AliEPSelectionTask::SetESDtrackCuts(TString status){ | |
356 | ||
357 | fStatus = status; | |
358 | if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); | |
359 | if (fStatus.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); | |
360 | fESDtrackCuts->SetPtRange(0.15,20); | |
361 | fESDtrackCuts->SetEtaRange(-0.8,0.8); | |
362 | } | |
363 | ||
364 | //__________________________________________________________________________ | |
08c34a44 | 365 | void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname) |
366 | { | |
ce7adfe9 | 367 | TFile f(infilename); |
368 | TObject* list = f.Get(listname); | |
369 | fPhiDist = (TH1F*)list->FindObject("fHOutPhi"); | |
370 | if (!fPhiDist) cout << "Phi Distribution not found!!!" << endl; | |
371 | ||
372 | Bool_t emptybins; | |
373 | ||
374 | int iter = 0; | |
375 | while (iter<3){ | |
376 | emptybins = kFALSE; | |
377 | ||
378 | for (int i=1; i<fPhiDist->GetNbinsX(); i++){ | |
379 | if (!((fPhiDist->GetBinContent(i))>0)) { | |
380 | emptybins = kTRUE; | |
381 | } | |
382 | } | |
383 | if (emptybins) { | |
384 | cout << "empty bins - rebinning!" << endl; | |
385 | fPhiDist->Rebin(); | |
386 | iter++; | |
387 | } | |
388 | else iter = 3; | |
389 | } | |
390 | ||
391 | if (emptybins) { | |
392 | AliError("After Maximum of rebinning still empty Phi-bins!!!"); | |
393 | } | |
394 | f.Close(); | |
395 | } | |
396 | ||
397 | //________________________________________________________________________ | |
398 | Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track) | |
399 | { | |
400 | Double_t ptweight=1; | |
401 | ||
402 | if (fUsePtWeight) { | |
403 | if (track->Pt()<2) ptweight=track->Pt(); | |
404 | else ptweight=2; | |
405 | } | |
406 | return ptweight*GetPhiWeight(track); | |
407 | } | |
408 | ||
409 | //________________________________________________________________________ | |
410 | Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track) | |
411 | { | |
412 | Double_t phiweight=1; | |
413 | ||
414 | if (fUsePhiWeight) { | |
415 | Double_t nParticles = fPhiDist->Integral(); | |
416 | Double_t nPhibins = fPhiDist->GetNbinsX(); | |
417 | ||
418 | Double_t Phi = track->Phi(); | |
419 | ||
420 | while (Phi<0) Phi += TMath::TwoPi(); | |
421 | while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi(); | |
422 | ||
08c34a44 | 423 | Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi())); |
ce7adfe9 | 424 | |
425 | if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue; | |
426 | } | |
427 | return phiweight; | |
428 | } |