]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/ESE/AliAnalysisTaskFemtoESE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / ESE / AliAnalysisTaskFemtoESE.cxx
CommitLineData
698fcf10 1#include <iostream>
2#include <math.h>
3#include "TChain.h"
4#include "TFile.h"
5#include "TKey.h"
6#include "TObject.h"
7#include "TObjString.h"
8#include "TList.h"
9#include "TTree.h"
10#include "TH1F.h"
11#include "TH1D.h"
12#include "TH2D.h"
13#include "TH3D.h"
14#include "THn.h"
15#include "THnSparse.h"
16#include "TProfile.h"
17#include "TProfile2D.h"
18#include "TCanvas.h"
19#include "TRandom3.h"
20#include "TF1.h"
21#include "TObjectTable.h"
22#include <vector>
23
e7199f47 24//#include "TStopwatch.h"
25
698fcf10 26#include "AliAnalysisTask.h"
27#include "AliAnalysisManager.h"
28
29
30#include "AliESDEvent.h"
31#include "AliESDInputHandler.h"
32#include "AliESDtrackCuts.h"
33
34#include "AliAODEvent.h"
35#include "AliAODInputHandler.h"
36#include "AliAODMCParticle.h"
37//#include "AliAnalysisUtils.h"
38#include "AliHelperPID.h"
39#include "AliEventPoolManager.h"
40
41#include "AliAnalysisTaskFemtoESE.h"
42
43//#include "AliSpectraAODEventCuts.h"
44//#include "AliSpectraAODTrackCuts.h"
45//#include "/opt/alice/aliroot/master/src/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.h"
46//#include "/opt/alice/aliroot/master/src/PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODTrackCuts.h"
47//#include "AliSpectraAODEventCuts.h"
48//#include "AliSpectraAODTrackCuts.h"
49
50ClassImp(AliAnalysisTaskFemtoESE)
51
52
53//________________________________________________________________________
54// Default constructor
55AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE() :
e7199f47 56AliAnalysisTaskSE(),
57 fAOD(0x0),
58 fOutputList(0x0),
59 fHelperPID(0x0),
60 fPoolMgr(0x0),
61 fEventCuts(0x0),
62 fTrackCuts(0x0),
63 fFilterBit(128),
64 fSelectBit(AliVEvent::kMB),
65 bIsLHC10h(1),
66 fEventCounter(0),
67 fMixingTracks(10000),
68 fBfield(0.),
69 fMinSepPairEta(0.),
70 fMinSepPairPhi(0.),
2d11620d 71 fQinvMin(-1.),
72 fMaxDcaXY(1000.),
73 fMaxDcaZ(1000.),
74 fPtMin(0.14),
75 fPtMax(1.5),
76 fEtaMax(0.8),
e7199f47 77 fShareQuality(0.5),
78 fShareFraction(0.05),
79 nCountSamePairs(0),
80 nCountMixedPairs(0),
81 nCountTracks(0),
82 fMinQPerc(-1000),
83 fMaxQPerc(1000),
d0570384 84 qlimit(0.2),
85 qbins(40),
e7199f47 86 fQPercDet(0),
87 fEPDet(0),
e7199f47 88 nKtBins(0),
89 nKtBins1(1),
90 ktBins(0),
91 nEPBins(0),
92 nEPBins1(1),
93 epBins(0),
f8caeea5 94 nEPBinsMix(0),
95 nEPBinsMix1(1),
96 epBinsMix(0),
e7199f47 97 nCentBins(0),
98 nCentBins1(1),
99 centBins(0),
100 nVzBins(0),
101 nVzBins1(1),
102 vzBins(0),
103 hq(0x0),
3b898cea 104 hqmix(0x0),
105 hqinv(0x0),
e7199f47 106 nqPercBinsLHC11h(1),
107 qPercBinsLHC11h(0),
108 hpx(0x0),
109 hpy(0x0),
110 hpz(0x0),
111 hpt(0x0),
112 hE(0x0),
113 hphieta(0x0),
114 hphieta_pid(0x0),
115 hpt_pid(0x0),
116 hvzcent(0x0),
117 hcent(0x0),
f8caeea5 118 hcentUnweighted(0x0),
e7199f47 119 hcentn(0x0),
120 hphistaretapair10(0x0),
121 hphistaretapair16(0x0),
122 hphistaretapair10a(0x0),
123 hphistaretapair16a(0x0),
124 hphistaretapair10b(0x0),
125 hphistaretapair16b(0x0),
126 hphietapair(0x0),
127 hphietapair2(0x0),
128 hpidid(0x0),
129 hkt(0x0),
130 hktcheck(0x0),
131 hkt3(0x0),
132 hdcaxy(0x0),
133 hdcaz(0x0),
134 hsharequal(0x0),
135 hsharefrac(0x0),
136 hsharequalmix(0x0),
137 hsharefracmix(0x0),
138 hPsiTPC(0x0),
139 hPsiV0A(0x0),
140 hPsiV0C(0x0),
f8caeea5 141 hShiftTPC(0x0),
142 hShiftV0A(0x0),
143 hShiftV0C(0x0),
144 hPsiMix(0x0),
e7199f47 145 hCheckEPA(0x0),
146 hCheckEPC(0x0),
147 hCheckEPmix(0x0),
f8caeea5 148 hAvDphi(0x0),
149 hNpairs(0x0),
150 hPairDphi(0x0),
151 hPairDphiMix(0x0),
e7199f47 152 hcentq(0x0),
f8caeea5 153 hMixedDistTracks(0x0),
154 hMixedDistEvents(0x0),
e7199f47 155 hQvecV0A(0x0),
156 hQvecV0C(0x0),
157 hresV0ATPC(0x0),
158 hresV0CTPC(0x0),
159 hresV0AV0C(0x0),
3b898cea 160 hqinvcheck(0x0),
e7199f47 161 hktbins(0x0),
162 hcentbins(0x0),
163 hepbins(0x0)
698fcf10 164{
165}
166
167//________________________________________________________________________
168// Constructor
169AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const char* name) :
170 AliAnalysisTaskSE(name),
171 fAOD(0x0),
172 fOutputList(0x0),
173 fHelperPID(0x0),
174 fPoolMgr(0x0),
175 fEventCuts(0x0),
176 fTrackCuts(0x0),
177 fFilterBit(128),
178 fSelectBit(AliVEvent::kMB),
76082c53 179 bIsLHC10h(1),
698fcf10 180 fEventCounter(0),
181 fMixingTracks(10000),
182 fBfield(0.),
183 fMinSepPairEta(0.),
184 fMinSepPairPhi(0.),
2d11620d 185 fQinvMin(-1.),
186 fMaxDcaXY(1000.),
187 fMaxDcaZ(1000.),
188 fPtMin(0.14),
189 fPtMax(1.5),
190 fEtaMax(0.8),
698fcf10 191 fShareQuality(0.5),
192 fShareFraction(0.05),
193 nCountSamePairs(0),
194 nCountMixedPairs(0),
195 nCountTracks(0),
196 fMinQPerc(-1000),
197 fMaxQPerc(1000),
d0570384 198 qlimit(0.2),
199 qbins(40),
698fcf10 200 fQPercDet(0),
201 fEPDet(0),
698fcf10 202 nKtBins(0),
203 nKtBins1(1),
204 ktBins(0),
205 nEPBins(0),
206 nEPBins1(1),
207 epBins(0),
f8caeea5 208 nEPBinsMix(0),
209 nEPBinsMix1(1),
210 epBinsMix(0),
698fcf10 211 nCentBins(0),
212 nCentBins1(1),
213 centBins(0),
214 nVzBins(0),
215 nVzBins1(1),
216 vzBins(0),
217 hq(0x0),
76082c53 218 hqmix(0x0),
3b898cea 219 hqinv(0x0),
76082c53 220 nqPercBinsLHC11h(1),
e7199f47 221 qPercBinsLHC11h(0),
222 hpx(0x0),
223 hpy(0x0),
224 hpz(0x0),
225 hpt(0x0),
226 hE(0x0),
227 hphieta(0x0),
228 hphieta_pid(0x0),
229 hpt_pid(0x0),
230 hvzcent(0x0),
231 hcent(0x0),
f8caeea5 232 hcentUnweighted(0x0),
e7199f47 233 hcentn(0x0),
234 hphistaretapair10(0x0),
235 hphistaretapair16(0x0),
236 hphistaretapair10a(0x0),
237 hphistaretapair16a(0x0),
238 hphistaretapair10b(0x0),
239 hphistaretapair16b(0x0),
240 hphietapair(0x0),
241 hphietapair2(0x0),
242 hpidid(0x0),
243 hkt(0x0),
244 hktcheck(0x0),
245 hkt3(0x0),
246 hdcaxy(0x0),
247 hdcaz(0x0),
248 hsharequal(0x0),
249 hsharefrac(0x0),
250 hsharequalmix(0x0),
251 hsharefracmix(0x0),
252 hPsiTPC(0x0),
253 hPsiV0A(0x0),
254 hPsiV0C(0x0),
f8caeea5 255 hShiftTPC(0x0),
256 hShiftV0A(0x0),
257 hShiftV0C(0x0),
258 hPsiMix(0x0),
e7199f47 259 hCheckEPA(0x0),
260 hCheckEPC(0x0),
261 hCheckEPmix(0x0),
f8caeea5 262 hAvDphi(0x0),
263 hNpairs(0x0),
264 hPairDphi(0x0),
265 hPairDphiMix(0x0),
e7199f47 266 hcentq(0x0),
f8caeea5 267 hMixedDistTracks(0x0),
268 hMixedDistEvents(0x0),
e7199f47 269 hQvecV0A(0x0),
270 hQvecV0C(0x0),
271 hresV0ATPC(0x0),
272 hresV0CTPC(0x0),
273 hresV0AV0C(0x0),
3b898cea 274 hqinvcheck(0x0),
e7199f47 275 hktbins(0x0),
276 hcentbins(0x0),
277 hepbins(0x0)
698fcf10 278{
279
280 Printf("*******************************************");
281 Printf("AliAnalysisTaskFemtoESE named %s",name);
282 Printf("*******************************************");
283
284 // default binning
aa5e3401 285 //SetEPBins(12,-TMath::Pi()/12.,2*TMath::Pi()-TMath::Pi()/12.);
d0570384 286 SetEPBins(12);
698fcf10 287 Double_t ktBinsTemp[5] = {0.2,0.3,0.4,0.5,0.7};
288 SetKtBins(4,ktBinsTemp);
aa5e3401 289 Double_t centBinsTemp[7] = {0,5,10,20,30,40,50};
290 SetCentBins(6,centBinsTemp);
291 Double_t vzBinsTemp[9] = {-8,-6,-4,-2,0,2,4,6,8};
292 SetVzBins(8,vzBinsTemp);
698fcf10 293
f8caeea5 294 nEPBinsMix = 6;
295 nEPBinsMix1 = 7;
296 epBinsMix = new Double_t[nEPBinsMix1];
297 for(Int_t ee = 0; ee < nEPBinsMix1; ee++)
298 {epBinsMix[ee] = (TMath::Pi()/(Double_t)nEPBinsMix)*((Double_t)ee)-TMath::Pi()/2.;}
299
698fcf10 300 vertex[0] = vertex[1] = vertex[2] = 0.;
301
302 DefineInput(0, TChain::Class());
303 DefineOutput(1, TList::Class());
304 DefineOutput(2, AliHelperPID::Class());
76082c53 305 DefineOutput(3, AliSpectraAODEventCuts::Class());
306 DefineOutput(4, AliSpectraAODTrackCuts::Class());
698fcf10 307
308}
309
310//________________________________________________________________________
311// destructor
312AliAnalysisTaskFemtoESE::~AliAnalysisTaskFemtoESE()
313{
314 /* do nothing yet */
315}
316//________________________________________________________________________
317// copy constructor
318AliAnalysisTaskFemtoESE::AliAnalysisTaskFemtoESE(const AliAnalysisTaskFemtoESE &/*obj*/) :
319 AliAnalysisTaskSE(),
320 fAOD(0x0),
321 fOutputList(0x0),
322 fHelperPID(0x0),
323 fPoolMgr(0x0),
324 fEventCuts(0x0),
325 fTrackCuts(0x0),
326 fFilterBit(128),
327 fSelectBit(AliVEvent::kMB),
76082c53 328 bIsLHC10h(1),
698fcf10 329 fEventCounter(0),
330 fMixingTracks(10000),
331 fBfield(0.),
332 fMinSepPairEta(0.),
333 fMinSepPairPhi(0.),
2d11620d 334 fQinvMin(-1.),
335 fMaxDcaXY(1000.),
336 fMaxDcaZ(1000.),
337 fPtMin(0.14),
338 fPtMax(1.5),
339 fEtaMax(0.8),
698fcf10 340 fShareQuality(0.5),
341 fShareFraction(0.05),
342 nCountSamePairs(0),
343 nCountMixedPairs(0),
344 nCountTracks(0),
345 fMinQPerc(-1000),
346 fMaxQPerc(1000),
d0570384 347 qlimit(0.2),
348 qbins(40),
698fcf10 349 fQPercDet(0),
350 fEPDet(0),
698fcf10 351 nKtBins(0),
352 nKtBins1(1),
353 ktBins(0),
354 nEPBins(0),
355 nEPBins1(1),
356 epBins(0),
f8caeea5 357 nEPBinsMix(0),
358 nEPBinsMix1(1),
359 epBinsMix(0),
698fcf10 360 nCentBins(0),
361 nCentBins1(1),
362 centBins(0),
363 nVzBins(0),
364 nVzBins1(1),
365 vzBins(0),
366 hq(0x0),
76082c53 367 hqmix(0x0),
3b898cea 368 hqinv(0x0),
76082c53 369 nqPercBinsLHC11h(1),
e7199f47 370 qPercBinsLHC11h(0),
371 hpx(0x0),
372 hpy(0x0),
373 hpz(0x0),
374 hpt(0x0),
375 hE(0x0),
376 hphieta(0x0),
377 hphieta_pid(0x0),
378 hpt_pid(0x0),
379 hvzcent(0x0),
380 hcent(0x0),
f8caeea5 381 hcentUnweighted(0x0),
e7199f47 382 hcentn(0x0),
383 hphistaretapair10(0x0),
384 hphistaretapair16(0x0),
385 hphistaretapair10a(0x0),
386 hphistaretapair16a(0x0),
387 hphistaretapair10b(0x0),
388 hphistaretapair16b(0x0),
389 hphietapair(0x0),
390 hphietapair2(0x0),
391 hpidid(0x0),
392 hkt(0x0),
393 hktcheck(0x0),
394 hkt3(0x0),
395 hdcaxy(0x0),
396 hdcaz(0x0),
397 hsharequal(0x0),
398 hsharefrac(0x0),
399 hsharequalmix(0x0),
400 hsharefracmix(0x0),
401 hPsiTPC(0x0),
402 hPsiV0A(0x0),
403 hPsiV0C(0x0),
f8caeea5 404 hShiftTPC(0x0),
405 hShiftV0A(0x0),
406 hShiftV0C(0x0),
407 hPsiMix(0x0),
e7199f47 408 hCheckEPA(0x0),
409 hCheckEPC(0x0),
410 hCheckEPmix(0x0),
f8caeea5 411 hAvDphi(0x0),
412 hNpairs(0x0),
413 hPairDphi(0x0),
414 hPairDphiMix(0x0),
e7199f47 415 hcentq(0x0),
f8caeea5 416 hMixedDistTracks(0x0),
417 hMixedDistEvents(0x0),
e7199f47 418 hQvecV0A(0x0),
419 hQvecV0C(0x0),
420 hresV0ATPC(0x0),
421 hresV0CTPC(0x0),
422 hresV0AV0C(0x0),
3b898cea 423 hqinvcheck(0x0),
e7199f47 424 hktbins(0x0),
425 hcentbins(0x0),
426 hepbins(0x0)
698fcf10 427{
428 /* do nothing yet */
429}
430//________________________________________________________________________
431// assignment operator
432AliAnalysisTaskFemtoESE& AliAnalysisTaskFemtoESE::operator=(const AliAnalysisTaskFemtoESE &/*obj*/)
433{
434 /* do nothing yet */
435 return *this;
436}
437
438//________________________________________________________________________
439void AliAnalysisTaskFemtoESE::UserCreateOutputObjects()
440{
441 // Create histograms
442
443 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
444 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
445 if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
e7199f47 446
698fcf10 447 fOutputList = new TList();
448 fOutputList->SetOwner();
449 fOutputList->SetName("fOutputList");
450
e7199f47 451 hpx = new TH1D("hpx","px",200,-2,2);
698fcf10 452 hpx->GetXaxis()->SetTitle("p_{x}");
453 fOutputList->Add(hpx);
e7199f47 454 hpy = new TH1D("hpy","py",200,-2,2);
698fcf10 455 hpy->GetXaxis()->SetTitle("p_{y}");
456 fOutputList->Add(hpy);
e7199f47 457 hpz = new TH1D("hpz","pz",200,-2,2);
698fcf10 458 hpz->GetXaxis()->SetTitle("p_{z}");
459 fOutputList->Add(hpz);
e7199f47 460 hpt = new TH1D("hpt","pt",100,0,2);
698fcf10 461 hpt->GetXaxis()->SetTitle("p_{t}");
462 fOutputList->Add(hpt);
e7199f47 463 hE = new TH1D("hE","E",100,0,2);
698fcf10 464 hE->GetXaxis()->SetTitle("E");
465 fOutputList->Add(hE);
e7199f47 466 hphieta = new TH2D("hphieta","track #varphi vs #eta",100,0,2*TMath::Pi(),80,-0.8,0.8);
698fcf10 467 hphieta->GetXaxis()->SetTitle("#varphi");
468 hphieta->GetYaxis()->SetTitle("#eta");
469 fOutputList->Add(hphieta);
e7199f47 470 hphieta_pid = new TH2D("hphieta_pid","PID check -- #Delta#varphi vs #Delta#eta",100,-0.3,0.3,100,-0.3,0.3);
698fcf10 471 hphieta_pid->GetXaxis()->SetTitle("#Delta#varphi");
472 hphieta_pid->GetYaxis()->SetTitle("#Delta#eta");
473 fOutputList->Add(hphieta_pid);
e7199f47 474 hpt_pid = new TH1D("hpt_pid","PID check -- #Delta p_{t}",100,-0.5,0.5);
698fcf10 475 hpt_pid->GetXaxis()->SetTitle("#Delta p_{t}");
476 fOutputList->Add(hpt_pid);
2d11620d 477 hvzcent = new TH2D("hvzcent","vz vs cent",20,-10,10,nCentBins,centBins);
698fcf10 478 hvzcent->GetXaxis()->SetTitle("v_{z}");
479 hvzcent->GetYaxis()->SetTitle("centrality");
480 fOutputList->Add(hvzcent);
3b898cea 481 hcent = new TH1D("hcent","cent",200,0,50);
698fcf10 482 hcent->GetXaxis()->SetTitle("centrality");
483 fOutputList->Add(hcent);
f8caeea5 484 hcentUnweighted = new TH1D("hcentUnweighted","cent - unweighted",200,0,50);
485 hcentUnweighted->GetXaxis()->SetTitle("centrality");
486 fOutputList->Add(hcentUnweighted);
e7199f47 487 hcentn = new TH2D("hcentn","cent vs npions",50,0,50,100,0,2000);
698fcf10 488 hcentn->GetXaxis()->SetTitle("Centrality");
489 hcentn->GetYaxis()->SetTitle("Number of pions");
490 fOutputList->Add(hcentn);
e7199f47 491 hphistaretapair10 = new TH3D("hphistaretapair10","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
698fcf10 492 hphistaretapair10->GetXaxis()->SetTitle("#Delta#varphi*");
493 hphistaretapair10->GetYaxis()->SetTitle("#Delta#eta");
494 hphistaretapair10->GetZaxis()->SetTitle("k_{T}");
495 fOutputList->Add(hphistaretapair10);
e7199f47 496 hphistaretapair16 = new TH3D("hphistaretapair16","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
698fcf10 497 hphistaretapair16->GetXaxis()->SetTitle("#Delta#varphi*");
498 hphistaretapair16->GetYaxis()->SetTitle("#Delta#eta");
499 hphistaretapair16->GetZaxis()->SetTitle("k_{T}");
500 fOutputList->Add(hphistaretapair16);
e7199f47 501 hphistaretapair10a = new TH3D("hphistaretapair10a","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
698fcf10 502 hphistaretapair10a->GetXaxis()->SetTitle("#Delta#varphi*");
503 hphistaretapair10a->GetYaxis()->SetTitle("#Delta#eta");
504 hphistaretapair10a->GetZaxis()->SetTitle("k_{T}");
505 fOutputList->Add(hphistaretapair10a);
e7199f47 506 hphistaretapair16a = new TH3D("hphistaretapair16a","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-0.15,0.15,100,-0.1,0.1,10,0,1);
698fcf10 507 hphistaretapair16a->GetXaxis()->SetTitle("#Delta#varphi*");
508 hphistaretapair16a->GetYaxis()->SetTitle("#Delta#eta");
509 hphistaretapair16a->GetZaxis()->SetTitle("k_{T}");
510 fOutputList->Add(hphistaretapair16a);
e7199f47 511 hphistaretapair10b = new TH3D("hphistaretapair10b","pair #Delta#varphi* vs #Delta#eta at r=1.0m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
698fcf10 512 hphistaretapair10b->GetXaxis()->SetTitle("#Delta#varphi*");
513 hphistaretapair10b->GetYaxis()->SetTitle("#Delta#eta");
514 hphistaretapair10b->GetZaxis()->SetTitle("k_{T}");
515 fOutputList->Add(hphistaretapair10b);
e7199f47 516 hphistaretapair16b = new TH3D("hphistaretapair16b","pair #Delta#varphi* vs #Delta#eta at r=1.6m",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
698fcf10 517 hphistaretapair16b->GetXaxis()->SetTitle("#Delta#varphi*");
518 hphistaretapair16b->GetYaxis()->SetTitle("#Delta#eta");
519 hphistaretapair16b->GetZaxis()->SetTitle("k_{T}");
520 fOutputList->Add(hphistaretapair16b);
e7199f47 521 hphietapair = new TH3D("hphietapair","pair #Delta#varphi vs #Delta#eta",100,-0.1,0.1,100,-0.1,0.1,10,0,1);
698fcf10 522 hphietapair->GetXaxis()->SetTitle("#Delta#varphi");
523 hphietapair->GetYaxis()->SetTitle("#Delta#eta");
524 hphietapair->GetZaxis()->SetTitle("k_{T}");
525 fOutputList->Add(hphietapair);
e7199f47 526 hphietapair2 = new TH3D("hphietapair2","pair #varphi vs #eta",100,-TMath::Pi(),TMath::Pi(),100,-1.6,1.6,10,0,1);
698fcf10 527 hphietapair2->GetXaxis()->SetTitle("#Delta#varphi");
528 hphietapair2->GetYaxis()->SetTitle("#eta");
529 hphietapair2->GetZaxis()->SetTitle("k_{T}");
530 fOutputList->Add(hphietapair2);
e7199f47 531 hpidid = new TH1D("hpidid","pid id",9,-4.5,4.5);
698fcf10 532 hpidid->GetXaxis()->SetTitle("track PID ID");
533 fOutputList->Add(hpidid);
f8caeea5 534 hkt = new TH2D("hkt","k_{T}",nCentBins,centBins,100,0,2);
535 hkt->GetXaxis()->SetTitle("centrality");
536 hkt->GetYaxis()->SetTitle("k_{T}");
698fcf10 537 fOutputList->Add(hkt);
e7199f47 538 hktcheck = new TH1D("hktcheck","k_{T} check",50,0,1);
698fcf10 539 hktcheck->GetXaxis()->SetTitle("k_{T}");
540 fOutputList->Add(hktcheck);
e7199f47 541 hkt3 = new TH3D("hkt3","kt vs pt",50,0,1,50,0,5,50,0,5);
698fcf10 542 hkt3->GetXaxis()->SetTitle("k_{T}");
543 hkt3->GetYaxis()->SetTitle("p_{T,1}");
544 hkt3->GetZaxis()->SetTitle("p_{T,2}");
545 fOutputList->Add(hkt3);
e7199f47 546 hdcaxy = new TH2D("hdcaxy","DCA xy",100,-5,5,100,-5,5);
698fcf10 547 hdcaxy->GetXaxis()->SetTitle("DCA x");
548 hdcaxy->GetYaxis()->SetTitle("DCA y");
549 fOutputList->Add(hdcaxy);
e7199f47 550 hdcaz = new TH1D("hdcaz","DCA z",100,-5,5);
698fcf10 551 hdcaz->GetXaxis()->SetTitle("DCA z");
552 fOutputList->Add(hdcaz);
e7199f47 553 hsharequal = new TH1D("hsharequal","Share Quality",102,-1.02,1.02);
698fcf10 554 hsharequal->GetXaxis()->SetTitle("Share Quality");
555 fOutputList->Add(hsharequal);
e7199f47 556 hsharefrac = new TH1D("hsharefrac","Share Fraction",100,0,1);
698fcf10 557 hsharefrac->GetXaxis()->SetTitle("Share Fraction");
558 fOutputList->Add(hsharefrac);
e7199f47 559 hsharequalmix = new TH1D("hsharequalmix","Share Quality -- mixed events",102,-1.02,1.02);
698fcf10 560 hsharequalmix->GetXaxis()->SetTitle("Share Quality");
561 fOutputList->Add(hsharequalmix);
e7199f47 562 hsharefracmix = new TH1D("hsharefracmix","Share Fraction -- mixed events",100,0,1);
698fcf10 563 hsharefracmix->GetXaxis()->SetTitle("Share Fraction");
564 fOutputList->Add(hsharefracmix);
f8caeea5 565 hPsiTPC = new TH2D("hPsiTPC","TPC EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
566 hPsiTPC->GetXaxis()->SetTitle("Centrality");
567 hPsiTPC->GetYaxis()->SetTitle("#Psi{TPC}");
698fcf10 568 fOutputList->Add(hPsiTPC);
f8caeea5 569 hPsiV0A = new TH2D("hPsiV0A","V0A EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
570 hPsiV0A->GetXaxis()->SetTitle("Centrality");
571 hPsiV0A->GetYaxis()->SetTitle("#Psi{V0A}");
698fcf10 572 fOutputList->Add(hPsiV0A);
f8caeea5 573 hPsiV0C = new TH2D("hPsiV0C","V0C EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
574 hPsiV0C->GetXaxis()->SetTitle("Centrality");
575 hPsiV0C->GetYaxis()->SetTitle("#Psi{V0C}");
698fcf10 576 fOutputList->Add(hPsiV0C);
f8caeea5 577 hShiftTPC = new TH2D("hShiftTPC","TPC shifting values",nCentBins,centBins,6,-0.5,5.5);
578 hShiftTPC->GetXaxis()->SetTitle("Centrality");
579 fOutputList->Add(hShiftTPC);
580 hShiftV0A = new TH2D("hShiftV0A","V0A shifting values",nCentBins,centBins,6,-0.5,5.5);
581 hShiftV0A->GetXaxis()->SetTitle("Centrality");
582 fOutputList->Add(hShiftV0A);
583 hShiftV0C = new TH2D("hShiftV0C","V0C shifting values",nCentBins,centBins,6,-0.5,5.5);
584 hShiftV0C->GetXaxis()->SetTitle("Centrality");
585 fOutputList->Add(hShiftV0C);
586 hPsiMix = new TH2D("hPsiMix","Mix EP",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
587 hPsiMix->GetXaxis()->SetTitle("Centrality");
588 hPsiMix->GetYaxis()->SetTitle("#Psi{Mix}");
589 fOutputList->Add(hPsiMix);
590 hCheckEPA = new TH2D("hCheckEPA","Check EP V0A",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
591 hCheckEPA->GetXaxis()->SetTitle("Centrality");
592 hCheckEPA->GetYaxis()->SetTitle("PsiV0A - PsiTPC");
698fcf10 593 fOutputList->Add(hCheckEPA);
f8caeea5 594 hCheckEPC = new TH2D("hCheckEPC","Check EP V0C",nCentBins,centBins,100,-1*TMath::Pi(),TMath::Pi());
595 hCheckEPC->GetXaxis()->SetTitle("Centrality");
596 hCheckEPC->GetYaxis()->SetTitle("PsiV0C - PsiTPC");
698fcf10 597 fOutputList->Add(hCheckEPC);
e7199f47 598 hCheckEPmix = new TH2D("hCheckEPmix","Check EP mixed events",100,-1*TMath::Pi(),TMath::Pi(),100,-1*TMath::Pi(),TMath::Pi());
698fcf10 599 hCheckEPmix->GetXaxis()->SetTitle("Psi1 - Psi_mix");
600 hCheckEPmix->GetYaxis()->SetTitle("Psi1 - Psi2");
601 fOutputList->Add(hCheckEPmix);
f8caeea5 602 hAvDphi = new TH3D("hAvDphi","average event plane angle",nEPBins,epBins,nCentBins,centBins,nKtBins,ktBins);
603 hAvDphi->GetXaxis()->SetTitle("EP bins");
604 hAvDphi->GetYaxis()->SetTitle("cent bins");
605 hAvDphi->GetZaxis()->SetTitle("kt bins");
606 fOutputList->Add(hAvDphi);
607 hNpairs = new TH3D("hNpairs","number of pairs",nEPBins,epBins,nCentBins,centBins,nKtBins,ktBins);
608 hNpairs->GetXaxis()->SetTitle("EP bins");
609 hNpairs->GetYaxis()->SetTitle("cent bins");
610 hNpairs->GetZaxis()->SetTitle("kt bins");
611 fOutputList->Add(hNpairs);
612 hPairDphi = new TH1D("hPairDphi","pairs wrt EP",100,0,2*TMath::Pi());
613 hPairDphi->GetXaxis()->SetTitle("#Delta#phi");
614 fOutputList->Add(hPairDphi);
615 hPairDphiMix = new TH1D("hPairDphiMix","mixed pairs wrt EP",100,0,2*TMath::Pi());
616 hPairDphiMix->GetXaxis()->SetTitle("#Delta#phi");
617 fOutputList->Add(hPairDphiMix);
e7199f47 618 hcentq = new TH2D("hcentq","qvec vs cent",100,0,100,50,0,50);
698fcf10 619 hcentq->GetXaxis()->SetTitle("q_{2} percentile");
620 hcentq->GetYaxis()->SetTitle("centrality");
621 fOutputList->Add(hcentq);
3b2ee55b 622 hMixedDistTracks = new TH2D("hMixedDistTracks", ";centrality;tracks", nCentBins, centBins, 200, 0, fMixingTracks * 2.);
f8caeea5 623 fOutputList->Add(hMixedDistTracks);
3b2ee55b 624 hMixedDistEvents = new TH2D("hMixedDistEvents", ";centrality;events", nCentBins, centBins, 41, -0.5, 40.5);
f8caeea5 625 fOutputList->Add(hMixedDistEvents);
e7199f47 626 hQvecV0A = new TH2D("hQvecV0A","Qvector in V0A",50,0,50,200,0,5);
76082c53 627 hQvecV0A->GetXaxis()->SetTitle("Centrality");
628 hQvecV0A->GetYaxis()->SetTitle("normalized Qvector");
629 fOutputList->Add(hQvecV0A);
e7199f47 630 hQvecV0C = new TH2D("hQvecV0C","Qvector in V0C",50,0,50,200,0,5);
76082c53 631 hQvecV0C->GetXaxis()->SetTitle("Centrality");
632 hQvecV0C->GetYaxis()->SetTitle("normalized Qvector");
633 fOutputList->Add(hQvecV0C);
698fcf10 634
aa5e3401 635 // resolution histograms
e7199f47 636 hresV0ATPC = new TH1D("hresV0ATPC","cent vs cos(2*(V0A-TPC))",nCentBins,centBins);
aa5e3401 637 hresV0ATPC->GetXaxis()->SetTitle("centrality");
698fcf10 638 fOutputList->Add(hresV0ATPC);
e7199f47 639 hresV0CTPC = new TH1D("hresV0CTPC","cent vs cos(2*(V0C-TPC))",nCentBins,centBins);
aa5e3401 640 hresV0CTPC->GetXaxis()->SetTitle("centrality");
698fcf10 641 fOutputList->Add(hresV0CTPC);
e7199f47 642 hresV0AV0C = new TH1D("hresV0AV0C","cent vs cos(2*(V0A-V0C))",nCentBins,centBins);
aa5e3401 643 hresV0AV0C->GetXaxis()->SetTitle("centrality");
698fcf10 644 fOutputList->Add(hresV0AV0C);
645
3b898cea 646 hqinvcheck = new TH3F("hqinvcheck","Qinv vs kt vs cent",100,0,1,50,0,1,10,0,50);
647 hqinvcheck->GetXaxis()->SetTitle("qinv");
648 hqinvcheck->GetYaxis()->SetTitle("kt");
649 hqinvcheck->GetZaxis()->SetTitle("centrality");
650 fOutputList->Add(hqinvcheck);
651
d0570384 652 Double_t limit1 = 8.71488;
653 Double_t limit2 = 12.0126;
3b2ee55b 654 hq = new TH3F****[2];
655 hqmix = new TH3F****[2];
656 hqinv = new TH3F****[2];
657 for(Int_t z = 0; z < 2; z++) // charge combinations
698fcf10 658 {
3b2ee55b 659 hq[z] = new TH3F***[nKtBins];
660 hqmix[z] = new TH3F***[nKtBins];
661 hqinv[z] = new TH3F***[nKtBins];
662 for(Int_t k = 0; k < nKtBins; k++)
698fcf10 663 {
3b2ee55b 664 hq[z][k] = new TH3F**[nEPBins];
665 hqmix[z][k] = new TH3F**[nEPBins];
666 hqinv[z][k] = new TH3F**[nEPBins];
667 for(Int_t e = 0; e < nEPBins; e++)
698fcf10 668 {
3b2ee55b 669 hq[z][k][e] = new TH3F*[nCentBins];
670 hqmix[z][k][e] = new TH3F*[nCentBins];
671 hqinv[z][k][e] = new TH3F*[nCentBins];
672 for(Int_t c = 0; c < nCentBins; c++)
673 {
d0570384 674 //hq[z][k][e][c] = new TH3F(Form("hq%i_k%i_e%i_c%i",z,k,e,c),Form("hq%i_k%i_e%i_c%i",z,k,e,c),60,-0.21,0.21,60,-0.21,0.21,60,-0.21,0.21);
675 hq[z][k][e][c] = new TH3F(Form("hq%i_k%i_e%i_c%i",z,k,e,c),Form("hq%i_k%i_e%i_c%i",z,k,e,c),qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit);
2d11620d 676 hq[z][k][e][c]->GetXaxis()->SetTitle("qout"); hq[z][k][e][c]->GetYaxis()->SetTitle("qside"); hq[z][k][e][c]->GetZaxis()->SetTitle("qlong");
d0570384 677 if(!(centBins[c]>limit2 || centBins[c+1]<limit1))
678 hq[z][k][e][c]->Sumw2();
679 // set sumw2 only for the correlation histograms which are filled with centrality weights (around cent=10)
3b2ee55b 680 fOutputList->Add(hq[z][k][e][c]);
d0570384 681 //hqmix[z][k][e][c] = new TH3F(Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),60,-0.21,0.21,60,-0.21,0.21,60,-0.21,0.21);
682 hqmix[z][k][e][c] = new TH3F(Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),Form("hqmix%i_k%i_e%i_c%i",z,k,e,c),qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit);
683 if(!(centBins[c]>limit2 || centBins[c+1]<limit1))
684 hqmix[z][k][e][c]->Sumw2();
3b2ee55b 685 fOutputList->Add(hqmix[z][k][e][c]);
d0570384 686 //hqinv[z][k][e][c] = new TH3F(Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),60,-0.21,0.21,60,-0.21,0.21,60,-0.21,0.21);
687 hqinv[z][k][e][c] = new TH3F(Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),Form("hqinv%i_k%i_e%i_c%i",z,k,e,c),qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit,qbins,-1.*qlimit,qlimit);
3b2ee55b 688 //hqinv[z][k][e][c]->Sumw2();
689 fOutputList->Add(hqinv[z][k][e][c]);
690 }
698fcf10 691 }
692 }
693 }
694
aa5e3401 695 // create dummy histograms which just hold the values of the kt, cent, ep bin edges
e7199f47 696 hktbins = new TH1F("hktbins","kt bins",nKtBins,ktBins);
698fcf10 697 fOutputList->Add(hktbins);
e7199f47 698 hcentbins = new TH1F("hcentbins","cent bins",nCentBins,centBins);
698fcf10 699 fOutputList->Add(hcentbins);
e7199f47 700 hepbins = new TH1F("hepbins","ep bins",nEPBins,epBins);
698fcf10 701 fOutputList->Add(hepbins);
698fcf10 702
698fcf10 703 Printf("************************");
f8caeea5 704 Printf("using the %s detector for event plane determination!",fEPDet ? "V0C" : "V0A");
705 Printf("using the %s detector for q-vector determination!",fQPercDet ? "V0C" : "V0A");
698fcf10 706 Printf("************************");
707
708 vertex[0] = vertex[1] = vertex[2] = 0.;
709
710 // event mixing pool
e7199f47 711 fPoolMgr = new AliEventPoolManager*[2];
698fcf10 712 Int_t poolsize = 1000;
f8caeea5 713 fPoolMgr[0] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
714 fPoolMgr[0]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
715 fPoolMgr[1] = new AliEventPoolManager(poolsize, fMixingTracks, nCentBins, centBins, nVzBins, vzBins,nEPBinsMix,epBinsMix);
716 fPoolMgr[1]->SetTargetValues(fMixingTracks, 0.01, 5); // check these values
698fcf10 717
718 nCountSamePairs = 0;
719 nCountMixedPairs = 0;
720 nCountTracks = 0;
721
e7199f47 722 //stopwatch = new TStopwatch();
723 //stopwatch->Start();
724
698fcf10 725 PostData(1, fOutputList);
726 PostData(2, fHelperPID);
76082c53 727 PostData(3, fEventCuts);
728 PostData(4, fTrackCuts);
e7199f47 729
698fcf10 730}
731
732//________________________________________________________________________
733void AliAnalysisTaskFemtoESE::UserExec(Option_t *)
e7199f47 734{
698fcf10 735 // Main loop
736 // Called for each event
737
738 //if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
739
740 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
741 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
742
743 //if(!EventCut(fAOD)) return;
744 if(!EventCut()) return;
745
746 fEventCounter++;
537c8bdb 747 if(fEventCounter%1000==0) Printf("=========== Event # %i ===========",fEventCounter);
698fcf10 748
749 AliCentrality *centrality;// for AODs and ESDs
750 const AliAODVertex *primaryVertexAOD;
751
752 //AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
753 //AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
754 //fPIDResponse = inputHandler->GetPIDResponse();
755
756 fBfield = fAOD->GetMagneticField();
757
758 /////////////////////////////////////////////////
759
760 Float_t centralityPercentile=0;
761
762 centrality = fAOD->GetCentrality();
763 centralityPercentile = centrality->GetCentralityPercentile("V0M");
764 if(centralityPercentile <= centBins[0]) return;
765 if(centralityPercentile > centBins[nCentBins]) return;
f8caeea5 766 Double_t centWeight = GetCentralityWeight(centralityPercentile);
698fcf10 767
768 primaryVertexAOD = fAOD->GetPrimaryVertex();
769 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
770 Double_t zvtx = vertex[2];
771 if(zvtx < vzBins[0] || zvtx > vzBins[nVzBins]) return; // Z-Vertex Cut
772 //cout<<"Centrality % = " << centralityPercentile << " z-vertex = " << zvtx << endl;
773
e7199f47 774 //stopwatch->Stop();
775 //Printf("%lf %lf",stopwatch->RealTime(),stopwatch->CpuTime());
776 //stopwatch->Start();
698fcf10 777
e7199f47 778 // get event plane from V0's
76082c53 779 if(!fEventCuts->IsSelected(fAOD,fTrackCuts)) {Printf("Error! Event not accepted by AliAODSpectraEventCuts!"); return;}
698fcf10 780 Double_t psiV0A = fEventCuts->GetPsiV0A();
781 Double_t psiV0C = fEventCuts->GetPsiV0C();
76082c53 782 Double_t qperc = -999;
f8caeea5 783 qperc = fEventCuts->GetQvecPercentile(fQPercDet);//0: VZERO-A 1: VZERO-C // now works for both LHC10h and LHC11h (for 0-50% centrality only!)
784 if(!bIsLHC10h)
76082c53 785 {
f8caeea5 786 //Printf("%lf %lf",qperc,GetQPercLHC11h(fEventCuts->GetqV0A()));
787 hQvecV0A->Fill(centralityPercentile,fEventCuts->GetqV0A(),centWeight);
788 hQvecV0C->Fill(centralityPercentile,fEventCuts->GetqV0C(),centWeight);
789 //Printf("q vector = %lf percentile = %lf",(fQPercDet==0 ? fEventCuts->GetqV0A() : fEventCuts->GetqV0C()),qperc);
76082c53 790 }
f8caeea5 791
698fcf10 792 if(psiV0A == -999) return;
793 if(psiV0C == -999) return;
794 if(qperc < fMinQPerc || qperc > fMaxQPerc) return;
795
f8caeea5 796 //ProcInfo_t procInfo;
797 //gSystem->GetProcInfo(&procInfo);
798 //Printf("beginning of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
799
e7199f47 800 Double_t psiEP = psiV0A;
801 if(fEPDet==1) psiEP = psiV0C;
698fcf10 802
803 Double_t sin2phi = 0, cos2phi = 0;
804
e7199f47 805 TObjArray* tracksPos = new TObjArray();
806 tracksPos->SetOwner(kTRUE);
807 TObjArray* tracksNeg = new TObjArray();
808 tracksNeg->SetOwner(kTRUE);
698fcf10 809
810 // Track loop -- select pions
811 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
f8caeea5 812 AliAODTrack* aodtrack = (AliAODTrack*)fAOD->GetTrack(i);
698fcf10 813 if (!aodtrack) continue;
814 if(!TrackCut(aodtrack)) continue;
815
d0570384 816 // check event plane angle using tracks in the TPC
817 if(aodtrack->Pt() < 2 && aodtrack->Pt() > 0.2)
818 {
819 sin2phi += (aodtrack->Pt())*sin(2*aodtrack->Phi());
820 cos2phi += (aodtrack->Pt())*cos(2*aodtrack->Phi());
821 }
822
698fcf10 823 // filter bit 7 PID method...
824 Int_t trackPID=999;
825 for(Int_t m = 0; m < fAOD->GetNumberOfTracks(); m++) {
f8caeea5 826 AliAODTrack* aodtrack2 = (AliAODTrack*)fAOD->GetTrack(m);
698fcf10 827 if (!aodtrack2) continue;
828 if(aodtrack->GetID() != (-aodtrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
829 trackPID=fHelperPID->GetParticleSpecies((AliVTrack*)aodtrack2,kTRUE);
830 hphieta_pid->Fill(aodtrack->Phi()-aodtrack2->Phi(),aodtrack->Eta()-aodtrack2->Eta());
831 hpt_pid->Fill(aodtrack->Pt()-aodtrack2->Pt());
832 //cout << aodtrack->Phi() << " " << aodtrack->Eta() << " " << aodtrack->Pt() << endl;
833 //cout << aodtrack2->Phi() << " " << aodtrack2->Eta() << " " << aodtrack2->Pt() << " " << dataID2 << endl;
834 }
835
836 hpidid->Fill((trackPID+1)*aodtrack->Charge());
837
838 // select pions
e7199f47 839 if(trackPID==0)
840 {
841 AliFemtoESEBasicParticle* particle = new AliFemtoESEBasicParticle(sqrt(pow(aodtrack->P(),2)+pow(0.13957, 2)),aodtrack->Px(),aodtrack->Py(),aodtrack->Pz(),aodtrack->Charge(),aodtrack->Phi(),aodtrack->Eta());
842 particle->SetPsiEP(psiEP);
843 particle->SetTPCClusterMap(aodtrack->GetTPCClusterMap());
844 particle->SetTPCSharedMap(aodtrack->GetTPCSharedMap());
845
846 if(particle->Charge()>0)
847 tracksPos->Add(particle);
848 if(particle->Charge()<0)
849 tracksNeg->Add(particle);
850
851 // track qa plots
852 hpx->Fill(particle->Px());
853 hpy->Fill(particle->Py());
854 hpz->Fill(particle->Pz());
855 hpt->Fill(particle->Pt());
856 hE->Fill(particle->E());
857 hphieta->Fill(particle->Phi(),particle->Eta());
858 }
698fcf10 859
860 }
861 // end track loop
862
e7199f47 863 Int_t ntracks = tracksPos->GetEntriesFast()+tracksNeg->GetEntriesFast();
698fcf10 864
865 // get EP from TPC, just to check
866 Double_t psiTPC = 0.;
d0570384 867 if(sin2phi != 0 && cos2phi != 0)
698fcf10 868 psiTPC = 0.5*atan2(sin2phi,cos2phi);
869 else return;
870
f8caeea5 871 hPsiTPC->Fill(centralityPercentile,psiTPC,centWeight);
872 hPsiV0A->Fill(centralityPercentile,psiV0A,centWeight);
873 hPsiV0C->Fill(centralityPercentile,psiV0C,centWeight);
698fcf10 874 Double_t dphiEP = psiTPC-psiV0A;
875 if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
876 if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
f8caeea5 877 hCheckEPA->Fill(centralityPercentile,dphiEP,centWeight);
698fcf10 878 dphiEP = psiTPC-psiV0C;
879 if(dphiEP>TMath::Pi()) dphiEP-=2*TMath::Pi();
880 if(dphiEP<-TMath::Pi()) dphiEP+=2*TMath::Pi();
f8caeea5 881 hCheckEPC->Fill(centralityPercentile,dphiEP,centWeight);
882
883 // values for EP shifting method
884 hShiftTPC->Fill(centralityPercentile,0.,cos(2*psiTPC)*centWeight);
885 hShiftTPC->Fill(centralityPercentile,1.,sin(2*psiTPC)*centWeight);
886 hShiftTPC->Fill(centralityPercentile,2.,cos(4*psiTPC)*centWeight);
887 hShiftTPC->Fill(centralityPercentile,3.,sin(4*psiTPC)*centWeight);
888 hShiftTPC->Fill(centralityPercentile,4.,cos(6*psiTPC)*centWeight);
889 hShiftTPC->Fill(centralityPercentile,5.,sin(6*psiTPC)*centWeight);
698fcf10 890
f8caeea5 891 hShiftV0A->Fill(centralityPercentile,0.,cos(2*psiV0A)*centWeight);
892 hShiftV0A->Fill(centralityPercentile,1.,sin(2*psiV0A)*centWeight);
893 hShiftV0A->Fill(centralityPercentile,2.,cos(4*psiV0A)*centWeight);
894 hShiftV0A->Fill(centralityPercentile,3.,sin(4*psiV0A)*centWeight);
895 hShiftV0A->Fill(centralityPercentile,4.,cos(6*psiV0A)*centWeight);
896 hShiftV0A->Fill(centralityPercentile,5.,sin(6*psiV0A)*centWeight);
698fcf10 897
f8caeea5 898 hShiftV0C->Fill(centralityPercentile,0.,cos(2*psiV0C)*centWeight);
899 hShiftV0C->Fill(centralityPercentile,1.,sin(2*psiV0C)*centWeight);
900 hShiftV0C->Fill(centralityPercentile,2.,cos(4*psiV0C)*centWeight);
901 hShiftV0C->Fill(centralityPercentile,3.,sin(4*psiV0C)*centWeight);
902 hShiftV0C->Fill(centralityPercentile,4.,cos(6*psiV0C)*centWeight);
903 hShiftV0C->Fill(centralityPercentile,5.,sin(6*psiV0C)*centWeight);
904
905
906 hcentq->Fill(qperc,centralityPercentile,centWeight);
698fcf10 907
698fcf10 908 nCountTracks += ntracks;
909 //cout << "Found " << ntracks << " pion tracks..." << endl;
910
f8caeea5 911 hvzcent->Fill(zvtx,centralityPercentile,centWeight);
912 hcentUnweighted->Fill(centralityPercentile);
913 hcent->Fill(centralityPercentile,centWeight);
914 hcentn->Fill(centralityPercentile,ntracks,centWeight);
698fcf10 915
916 // resolution histograms
f8caeea5 917 hresV0ATPC->Fill(centralityPercentile,cos(2*(psiV0A-psiTPC))*centWeight);
918 hresV0CTPC->Fill(centralityPercentile,cos(2*(psiV0C-psiTPC))*centWeight);
919 hresV0AV0C->Fill(centralityPercentile,cos(2*(psiV0A-psiV0C))*centWeight);
698fcf10 920
f8caeea5 921 AliEventPool* poolPos = fPoolMgr[0]->GetEventPool(centralityPercentile,zvtx,psiEP);
922 if (!poolPos) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
923 //Printf("Positive pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolPos->GetCurrentNEvents());
924 AliEventPool* poolNeg = fPoolMgr[1]->GetEventPool(centralityPercentile,zvtx,psiEP);
925 if (!poolNeg) AliFatal(Form("No pool found for centrality = %f, vz = %f, ep = %f", centralityPercentile, zvtx,psiEP));
926 //Printf("Negative pool found for centrality = %f, vz = %f, ep = %f with %i events in it", centralityPercentile, zvtx,psiEP,poolNeg->GetCurrentNEvents());
698fcf10 927
3b2ee55b 928 TrackLoop(tracksPos,poolPos,0,psiEP,centralityPercentile);
929 TrackLoop(tracksNeg,poolNeg,1,psiEP,centralityPercentile);
e7199f47 930 //TObjArray* clonedtracks = CloneAndReduceTrackList(tracks,psiEP);
931 poolPos->UpdatePool(tracksPos);
932 poolNeg->UpdatePool(tracksNeg);
933 //cout << "pool contains " << pool->GetCurrentNEvents() << " events and " << pool->NTracksInPool() << " tracks." << endl;
934 //tracks->Clear();
935
936 //delete tracks;
937
f8caeea5 938 //gSystem->GetProcInfo(&procInfo);
939 //Printf("end of event: ResMem %ld VMem %ld", procInfo.fMemResident, procInfo.fMemVirtual);
940
e7199f47 941 // Post output data.
942 PostData(1, fOutputList);
943 PostData(2, fHelperPID);
944 PostData(3, fEventCuts);
945 PostData(4, fTrackCuts);
946}
947
948
3b2ee55b 949void AliAnalysisTaskFemtoESE::TrackLoop(TObjArray *tracks, AliEventPool *pool, Int_t z, Double_t psiEP, Float_t centralityPercentile)
e7199f47 950{
3b2ee55b 951 // z=0 for positive charges, z=1 for negative charges
e7199f47 952 Double_t kt = 0;
953 Double_t qout=0, qside=0, qlong=0;
954 Double_t pVect1[4] = {0,0,0,0};
955 Double_t pVect2[4] = {0,0,0,0};
956 Int_t k, e, c; //bin indices for histograms
957
f8caeea5 958 Double_t centWeight = GetCentralityWeight(centralityPercentile);
959
e7199f47 960 Int_t ntracks = tracks->GetEntriesFast();
698fcf10 961
f8caeea5 962 Int_t countMixedEvents = 0, countMixedTracks = 0;
963
964 // same event loop
698fcf10 965 for(Int_t j = 0; j < ntracks; j++)
966 {
967 //cout << endl << j << " ";
e7199f47 968 AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
969 pVect1[0]=track1->E();
698fcf10 970 pVect1[1]=track1->Px();
971 pVect1[2]=track1->Py();
972 pVect1[3]=track1->Pz();
973 //cout << pVect1[0] << " " << pVect1[1] << " " << pVect1[2] << " " << pVect1[3] << endl;
974
698fcf10 975 for(Int_t i = j+1; i < ntracks; i++)
976 {
e7199f47 977 AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)tracks->At(i);
978
979 Double_t deltaphistar10 = DeltaPhiStar(track1,track2,1.0);
980 Double_t deltaphistar16 = DeltaPhiStar(track1,track2,1.6);
698fcf10 981
e7199f47 982 hphistaretapair10->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
983 hphistaretapair16->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
698fcf10 984
985 if(!PairCut(track1,track2,kFALSE)) continue;
986
e7199f47 987 hphistaretapair10a->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
988 hphistaretapair16a->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
989 hphistaretapair10b->Fill(deltaphistar10,track1->Eta()-track2->Eta(),kt);
990 hphistaretapair16b->Fill(deltaphistar16,track1->Eta()-track2->Eta(),kt);
698fcf10 991
e7199f47 992 pVect2[0]=track2->E();
698fcf10 993 pVect2[1]=track2->Px();
994 pVect2[2]=track2->Py();
995 pVect2[3]=track2->Pz();
996
997 //qinv = GetQinv(pVect1, pVect2); // = qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
998 GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
d0570384 999 if(!(fabs(qout)<qlimit && fabs(qside)<qlimit && fabs(qlong)<qlimit)) continue;
1000 nCountSamePairs++;
698fcf10 1001 kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
f8caeea5 1002 hkt->Fill(centralityPercentile,kt,centWeight);
698fcf10 1003 hkt3->Fill(kt,track1->Pt(),track2->Pt());
1004 Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEP); // angle to event plane in correct range
aa5e3401 1005 if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
698fcf10 1006 hktcheck->Fill(kt);
d0570384 1007 if(deltaphi > TMath::Pi())
1008 {
1009 hq[z][k][e][c]->Fill(qout,-1.*qside,-1.*qlong,centWeight);
1010 hqinv[z][k][e][c]->Fill(qout,-1.*qside,-1.*qlong,GetQinv(pVect1, pVect2)*centWeight);
1011 }
1012 else
1013 {
1014 hq[z][k][e][c]->Fill(qout,qside,qlong,centWeight);
1015 hqinv[z][k][e][c]->Fill(qout,qside,qlong,GetQinv(pVect1, pVect2)*centWeight);
1016 }
f8caeea5 1017 hNpairs->Fill(deltaphi,centralityPercentile,kt);
1018 hAvDphi->Fill(deltaphi,centralityPercentile,kt,deltaphi);
1019 hPairDphi->Fill(deltaphi);
1020 hqinvcheck->Fill(GetQinv(pVect1, pVect2),kt,centralityPercentile,centWeight);
698fcf10 1021 Double_t dphi = track1->Phi()-track2->Phi();
1022 if(dphi<-TMath::Pi()) dphi += 2*TMath::Pi();
1023 if(dphi>TMath::Pi()) dphi -= 2*TMath::Pi();
1024 hphietapair->Fill(dphi,track1->Eta()-track2->Eta(),kt);
1025 hphietapair2->Fill(dphi,track1->Eta()-track2->Eta(),kt);
1026 //cout << k << " ";
1027 }
f8caeea5 1028 }
698fcf10 1029
f8caeea5 1030 // mixed event loop
1031 countMixedTracks = 0;
1032 countMixedEvents = 0;
1033 if (pool->IsReady())
1034 {
1035 Int_t nMix = pool->GetCurrentNEvents();
1036
1037 for (Int_t jMix=0; jMix<nMix; jMix++) // loop over events in pool
698fcf10 1038 {
f8caeea5 1039 TObjArray* bgTracks = pool->GetEvent(jMix);
1040 Int_t ntracksmix = bgTracks->GetEntriesFast();
1041 if(ntracksmix <= 0) continue;
1042
1043 // compare event planes of two events
1044 AliFemtoESEBasicParticle* tracktest = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(0);
1045 Double_t psiEPmixtemp = tracktest->GetPsiEP();
1046 /*Double_t dphiEPtest = fPsiEPmixtemp-psiEP;
1047 while(dphiEPtest>2*TMath::Pi()) dphiEPtest-=2*TMath::Pi();
1048 while(dphiEPtest<0) dphiEPtest+=2*TMath::Pi();
1049 if(dphiEPtest>TMath::Pi()) dphiEPtest-=TMath::Pi();
1050 if(dphiEPtest>TMath::Pi()/2.) dphiEPtest = TMath::Pi()-dphiEPtest;
1051 if(dphiEPtest > TMath::Pi()/6.) continue; // event planes must be within pi/6
1052 */
1053
1054 Double_t psiEPmix = 0.5*atan2(sin(2*psiEP)+sin(2*psiEPmixtemp),cos(2*psiEP)+cos(2*psiEPmixtemp));
1055 hPsiMix->Fill(centralityPercentile,psiEPmix,centWeight);
1056
1057 Double_t dphimix = psiEP-psiEPmix;
1058 if(dphimix < -TMath::Pi()) dphimix += 2*TMath::Pi();
1059 if(dphimix > TMath::Pi()) dphimix -= 2*TMath::Pi();
1060 Double_t dphi12 = psiEP-psiEPmixtemp;
1061 if(dphi12 < -TMath::Pi()) dphi12 += 2*TMath::Pi();
1062 if(dphi12 > TMath::Pi()) dphi12 -= 2*TMath::Pi();
1063 hCheckEPmix->Fill(dphimix,dphi12);
1064
1065 countMixedTracks += ntracksmix;
1066 countMixedEvents += 1;
1067
1068 for(Int_t j = 0; j < ntracks; j++)
698fcf10 1069 {
f8caeea5 1070 //cout << endl << j << " ";
1071 AliFemtoESEBasicParticle* track1 = (AliFemtoESEBasicParticle*)tracks->At(j);
1072 pVect1[0]=track1->E();
1073 pVect1[1]=track1->Px();
1074 pVect1[2]=track1->Py();
1075 pVect1[3]=track1->Pz();
1076
698fcf10 1077 for(Int_t i = 0; i < ntracksmix; i++)
1078 {
1079 AliFemtoESEBasicParticle* track2 = (AliFemtoESEBasicParticle*)bgTracks->UncheckedAt(i);
1080
1081 if(!PairCut(track1,track2,kTRUE)) continue;
1082
1083 pVect2[0]=track2->E();
1084 pVect2[1]=track2->Px();
1085 pVect2[2]=track2->Py();
1086 pVect2[3]=track2->Pz();
1087
f8caeea5 1088 if(psiEPmixtemp != track2->GetPsiEP()) AliFatal("Error! Event plane angles are wrong in mixing!!");
698fcf10 1089
1090 //qinv = GetQinv(pVect1, pVect2); // qinv**2 = (P1x-P2x)**2 + (P1y-P2y)**2 + (P1z-P2z)**2 - (P1t-P2t)**2
1091 GetQosl(pVect1, pVect2, qout, qside, qlong); // qout, qside, qlong = components of Q=P1-P2 in the P=P1+P2 frame
d0570384 1092 if(!(fabs(qout)<qlimit && fabs(qside)<qlimit && fabs(qlong)<qlimit)) continue;
1093 nCountMixedPairs++;
698fcf10 1094 kt = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.; // = Kt = |pT1+pT2|/2
f8caeea5 1095 Double_t deltaphi = GetDeltaPhiEP(pVect1[1],pVect1[2],pVect2[1],pVect2[2],psiEPmix); // angle to event plane in correct range
698fcf10 1096 //Double_t weight = 1./(Double_t)nMix;
aa5e3401 1097 if(!FindBin(kt,deltaphi,centralityPercentile,k,e,c)) continue;
d0570384 1098 if(deltaphi > TMath::Pi())
1099 hqmix[z][k][e][c]->Fill(qout,-1.*qside,-1.*qlong,centWeight);
1100 else
1101 hqmix[z][k][e][c]->Fill(qout,qside,qlong,centWeight);
f8caeea5 1102
1103 hPairDphiMix->Fill(deltaphi);
698fcf10 1104
1105 }
1106 }
698fcf10 1107 }
1108 }
1109
f8caeea5 1110 hMixedDistTracks->Fill(centralityPercentile, countMixedTracks);
1111 hMixedDistEvents->Fill(centralityPercentile, countMixedEvents);
1112 //Printf("mixed tracks: %i mixed events: %i",countMixedTracks,countMixedEvents);
698fcf10 1113}
e7199f47 1114
1115
1116
1117
1118
698fcf10 1119//________________________________________________________________________
1120void AliAnalysisTaskFemtoESE::Terminate(Option_t *)
1121{
1122
1123 if(ktBins) delete [] ktBins;
1124 if(epBins) delete [] epBins;
1125 if(centBins) delete [] centBins;
1126 if(vzBins) delete [] vzBins;
76082c53 1127 if(qPercBinsLHC11h) delete [] qPercBinsLHC11h;
698fcf10 1128
1129 // Called once at the end of the query
1130
537c8bdb 1131 Printf("Done");
698fcf10 1132
1133}
1134
1135
1136/*
1137
1138//________________________________________________________________________
1139Bool_t AliAnalysisTaskFemtoESE::AcceptPair(AliChaoticityTrackStruct *first, AliChaoticityTrackStruct *second)
1140{
1141
e7199f47 1142if(fabs(first->fEta-second->fEta) > fMinSepPairEta) return kTRUE;
698fcf10 1143
e7199f47 1144// propagate through B field to r=1m
1145Float_t phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.15/first->fPt);// 0.15 for D=1m
1146if(phi1 > 2*PI) phi1 -= 2*PI;
1147if(phi1 < 0) phi1 += 2*PI;
1148Float_t phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.15/second->fPt);// 0.15 for D=1m
1149if(phi2 > 2*PI) phi2 -= 2*PI;
1150if(phi2 < 0) phi2 += 2*PI;
698fcf10 1151
e7199f47 1152Float_t deltaphi = phi1 - phi2;
1153if(deltaphi > PI) deltaphi -= 2*PI;
1154if(deltaphi < -PI) deltaphi += 2*PI;
1155deltaphi = fabs(deltaphi);
698fcf10 1156
e7199f47 1157if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
698fcf10 1158
1159
e7199f47 1160// propagate through B field to r=1.6m
1161phi1 = first->fPhi - asin(first->fCharge*(0.1*fBfield)*0.24/first->fPt);// mine. 0.24 for D=1.6m
1162if(phi1 > 2*PI) phi1 -= 2*PI;
1163if(phi1 < 0) phi1 += 2*PI;
1164phi2 = second->fPhi - asin(second->fCharge*(0.1*fBfield)*0.24/second->fPt);// mine. 0.24 for D=1.6m
1165if(phi2 > 2*PI) phi2 -= 2*PI;
1166if(phi2 < 0) phi2 += 2*PI;
698fcf10 1167
e7199f47 1168deltaphi = phi1 - phi2;
1169if(deltaphi > PI) deltaphi -= 2*PI;
1170if(deltaphi < -PI) deltaphi += 2*PI;
1171deltaphi = fabs(deltaphi);
698fcf10 1172
e7199f47 1173if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
698fcf10 1174
1175
1176
e7199f47 1177//
698fcf10 1178
e7199f47 1179Int_t ncl1 = first->fClusterMap.GetNbits();
1180Int_t ncl2 = second->fClusterMap.GetNbits();
1181Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1182Double_t shfrac = 0; Double_t qfactor = 0;
1183for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
1184if (first->fClusterMap.TestBitNumber(imap) && second->fClusterMap.TestBitNumber(imap)) {// Both clusters
1185if (first->fSharedMap.TestBitNumber(imap) && second->fSharedMap.TestBitNumber(imap)) { // Shared
1186sumQ++;
1187sumCls+=2;
1188sumSha+=2;}
1189else {sumQ--; sumCls+=2;}
1190}
1191else if (first->fClusterMap.TestBitNumber(imap) || second->fClusterMap.TestBitNumber(imap)) {// Non shared
1192sumQ++;
1193sumCls++;}
1194}
1195if (sumCls>0) {
1196qfactor = sumQ*1.0/sumCls;
1197shfrac = sumSha*1.0/sumCls;
1198}
698fcf10 1199
e7199f47 1200if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
698fcf10 1201
1202
e7199f47 1203return kTRUE;
698fcf10 1204
1205
1206}
1207//________________________________________________________________________
1208Float_t AliAnalysisTaskFemtoESE::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
1209{
e7199f47 1210Float_t arg = G_Coeff/qinv;
698fcf10 1211
e7199f47 1212if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
1213else {return (exp(-arg)-1)/(-arg);}
698fcf10 1214
1215}
1216//________________________________________________________________________
1217void AliAnalysisTaskFemtoESE::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
1218{
e7199f47 1219Int_t j, k;
1220Int_t a = i2 - i1;
1221for (Int_t i = i1; i < i2+1; i++) {
1222j = (Int_t) (gRandom->Rndm() * a);
1223k = iarr[j];
1224iarr[j] = iarr[i];
1225iarr[i] = k;
1226}
698fcf10 1227}
1228*/
1229//________________________________________________________________________
f8caeea5 1230Double_t AliAnalysisTaskFemtoESE::GetQinv(Double_t track1[], Double_t track2[]){
698fcf10 1231
e7199f47 1232 Double_t qinv2=1.0;
f8caeea5 1233 qinv2 = pow(track1[1]-track2[1],2.) + pow(track1[2]-track2[2],2.) + pow(track1[3]-track2[3],2.) - pow(track1[0]-track2[0],2.);
1234 if(qinv2 >= 0.) return sqrt(qinv2);
1235 else return -1.*sqrt(-1.*qinv2);
698fcf10 1236}
1237
1238//________________________________________________________________________
1239void AliAnalysisTaskFemtoESE::GetQosl(Double_t track1[], Double_t track2[], Double_t& qout, Double_t& qside, Double_t& qlong){
1240
1241 Double_t p0 = track1[0] + track2[0];
1242 Double_t px = track1[1] + track2[1];
1243 Double_t py = track1[2] + track2[2];
1244 Double_t pz = track1[3] + track2[3];
1245
1246 Double_t mt = sqrt(p0*p0 - pz*pz);
1247 Double_t pt = sqrt(px*px + py*py);
1248
1249 Double_t v0 = track1[0] - track2[0];
1250 Double_t vx = track1[1] - track2[1];
1251 Double_t vy = track1[2] - track2[2];
1252 Double_t vz = track1[3] - track2[3];
1253
1254 if(gRandom->Rndm()<0.5)
1255 {
1256 v0 = -v0;
1257 vx = -vx;
1258 vy = -vy;
1259 vz = -vz;
1260 }
1261
1262 //cout << p0 << " " << px << " " << py << " " << pz << " " << v0 << " " << vx << " " << vy << " " << vz << " " << mt << " " << pt << endl;
1263
1264 qout = (px*vx + py*vy)/pt;
1265 qside = (px*vy - py*vx)/pt;
1266 qlong = (p0*vz - pz*v0)/mt;
1267}
1268
1269//________________________________________________________________________
1270Bool_t AliAnalysisTaskFemtoESE::EventCut(/*AliAODEvent* fevent*/){
1271
1272 // Trigger Cut
1273
1274 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1275 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1276 Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1277 if(!isSelected1 && !isSelected2 && !isSelected3) {return kFALSE;}
1278
e7199f47 1279 /*
1280 // Pile-up rejection
1281 AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1282 if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1283 else AnaUtil->SetUseMVPlpSelection(kFALSE);
1284 Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1285 if(pileUpCase) return;
1286
1287 // Vertexing
1288 primaryVertexAOD = fAOD->GetPrimaryVertex();
1289 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
698fcf10 1290
e7199f47 1291 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1292 ((TH3D*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
698fcf10 1293
e7199f47 1294 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1295 AliAODTrack* aodtrack = fAOD->GetTrack(i);
1296 if (!aodtrack) continue;
1297 AODTracks++;
1298 if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1299 FBTracks++;
1300 }
1301 ((TH1D*)fOutputList->FindObject("fMultDist2"))->Fill(FBTracks);
698fcf10 1302
e7199f47 1303 //if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Old Pile-up cut
1304 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
698fcf10 1305
e7199f47 1306 ((TH1D*)fOutputList->FindObject("fMultDist3"))->Fill(FBTracks);
698fcf10 1307
e7199f47 1308 fBfield = fAOD->GetMagneticField();
698fcf10 1309
e7199f47 1310 for(Int_t i=0; i<fZvertexBins; i++){
1311 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1312 zbin=i;
1313 break;
1314 }
1315 }
1316 */
698fcf10 1317
e7199f47 1318 return kTRUE;
698fcf10 1319
1320}
1321
1322//________________________________________________________________________
1323Bool_t AliAnalysisTaskFemtoESE::TrackCut(AliAODTrack* ftrack){
1324
1325 if (!ftrack->TestFilterBit(fFilterBit)) return kFALSE;
1326
2d11620d 1327 if(ftrack->Pt() < fPtMin) return kFALSE;
1328 if(ftrack->Pt() > fPtMax) return kFALSE;
1329 if(fabs(ftrack->Eta()) > fEtaMax) return kFALSE;
698fcf10 1330
1331 if(ftrack->GetTPCNcls() < 80) return kFALSE;// TPC nCluster cut
1332
2d11620d 1333 //Double_t trackdca[3] = {ftrack->XAtDCA(),ftrack->YAtDCA(),ftrack->ZAtDCA()};
1334 //Double_t dcaxy = sqrt( pow(trackdca[0] - vertex[0],2) + pow(trackdca[1] - vertex[1],2));
1335 //Double_t dcaz = sqrt( pow(trackdca[2] - vertex[2],2));
1336 Double_t dcaxy = sqrt( pow(ftrack->XAtDCA(),2) + pow(ftrack->YAtDCA(),2));
1337 Double_t dcaz = fabs(ftrack->ZAtDCA());
1338 if(dcaxy > fMaxDcaXY) return kFALSE;
1339 if(dcaz > fMaxDcaZ) return kFALSE;
1340 hdcaxy->Fill(ftrack->XAtDCA(),ftrack->YAtDCA());
1341 hdcaz->Fill(ftrack->ZAtDCA());
698fcf10 1342
1343
e7199f47 1344 //// FilterBit Overlap Check
1345 //if(fFilterBit != 7){
1346 // Bool_t goodTrackOtherFB = kFALSE;
1347 // if(fMCcase && fAOD->GetRunNumber()<=126437) goodTrackOtherFB=kTRUE;// FB7 to FB5 mapping in 10f6a MC does not work
1348 //
1349 // for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1350 // AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1351 // if(!aodtrack2) continue;
1352 // if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1353 //
1354 // if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1355 //
1356 // }
1357 // if(!goodTrackOtherFB) continue;
698fcf10 1358 // }
1359
1360 return kTRUE;
1361
1362}
1363
1364//________________________________________________________________________
e7199f47 1365Bool_t AliAnalysisTaskFemtoESE::PairCut(AliFemtoESEBasicParticle* ftrack1, AliFemtoESEBasicParticle* ftrack2, Bool_t mix){
698fcf10 1366
1367 // check for same charge
1368 if(ftrack2->Charge() != ftrack1->Charge()) return kFALSE;
1369
1370 // qinv cut
2d11620d 1371 if(fQinvMin > 0)
1372 {
1373 Double_t trackvec1[4] = {ftrack1->E(),ftrack1->Px(),ftrack1->Py(),ftrack1->Pz()};
1374 Double_t trackvec2[4] = {ftrack2->E(),ftrack2->Px(),ftrack2->Py(),ftrack2->Pz()};
1375 Double_t qinv = GetQinv(trackvec1,trackvec2);
1376 if(qinv < fQinvMin) return kFALSE; // qinv < 0.005
1377 }
698fcf10 1378 // deltaEta x deltaPhi* cut
1379 if(fabs(ftrack1->Eta()-ftrack2->Eta()) < fMinSepPairEta)
1380 {
2d11620d 1381 //Double_t deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.0); // angular separation at r=1m
1382 //deltaphistar = fabs(deltaphistar);
1383 //if(deltaphistar < fMinSepPairPhi) return kFALSE;
1384 Double_t deltaphistar = DeltaPhiStar(ftrack1,ftrack2,1.6); // angular separation at r=1.6m
698fcf10 1385 deltaphistar = fabs(deltaphistar);
1386 if(deltaphistar < fMinSepPairPhi) return kFALSE;
1387 }
1388
1389 // share fraction & share quality cut
1390 TBits clusterMap1 = (TBits)(ftrack1->GetTPCClusterMap());
1391 TBits sharedMap1 = (TBits)(ftrack1->GetTPCSharedMap());
1392 TBits clusterMap2 = (TBits)(ftrack2->GetTPCClusterMap());
1393 TBits sharedMap2 = (TBits)(ftrack2->GetTPCSharedMap());
1394
1395 Int_t ncl1 = clusterMap1.GetNbits();
1396 Int_t ncl2 = clusterMap2.GetNbits();
1397 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
1398 Double_t shfrac = 0; Double_t qfactor = 0;
1399 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++)
1400 {
1401 if (clusterMap1.TestBitNumber(imap) && clusterMap2.TestBitNumber(imap)) // Both clusters
1402 {
1403 if (sharedMap1.TestBitNumber(imap) && sharedMap2.TestBitNumber(imap)) // Shared
1404 {
1405 sumQ++;
1406 sumCls+=2;
1407 sumSha+=2;}
1408 else {sumQ--; sumCls+=2;}
1409 }
1410 else if (clusterMap1.TestBitNumber(imap) || clusterMap2.TestBitNumber(imap)) // Non shared
1411 {
1412 sumQ++;
1413 sumCls++;
1414 }
1415 }
1416 if (sumCls>0)
1417 {
1418 qfactor = sumQ*1.0/sumCls;
1419 shfrac = sumSha*1.0/sumCls;
1420 }
1421
1422 // sumCls -- number of clusters in track 1 + number of clusters in track 2 (clusters in both tracks counted twice)
1423 // sumSha -- number of shared clusters (counted twice)
1424 // sumQ -- ?
698fcf10 1425
1426 if(!mix)
1427 {
e7199f47 1428 hsharequal->Fill(qfactor);
1429 hsharefrac->Fill(shfrac);
698fcf10 1430 }
1431 else
1432 {
e7199f47 1433 hsharequalmix->Fill(qfactor);
1434 hsharefracmix->Fill(shfrac);
698fcf10 1435 }
e7199f47 1436
698fcf10 1437
1438 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
1439
1440 return kTRUE;
1441}
1442
1443Double_t AliAnalysisTaskFemtoESE::DeltaPhiStar(AliAODTrack* ftrack1, AliAODTrack* ftrack2, Double_t r)
1444{
1445 Double_t phi1 = ftrack1->Phi() - asin(0.3*ftrack1->Charge()*(0.1*fBfield)*r/(2.*ftrack1->Pt())); // magnetic field converted from kilogauss to tesla
1446 Double_t phi2 = ftrack2->Phi() - asin(0.3*ftrack2->Charge()*(0.1*fBfield)*r/(2.*ftrack2->Pt()));
1447
1448 Double_t deltaphi = phi1 - phi2;
1449 while(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1450 while(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1451
1452 return deltaphi;
1453}
1454
e7199f47 1455/*TObjArray* AliAnalysisTaskFemtoESE::CloneAndReduceTrackList(TObjArray* tracks, Double_t psi)
1456 {
698fcf10 1457 // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1458
1459 TObjArray* tracksClone = new TObjArray;
1460 tracksClone->SetOwner(kTRUE);
1461
1462 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1463 {
e7199f47 1464 AliAODTrack* particle = (AliAODTrack*) tracks->UncheckedAt(i);
1465 AliFemtoESEBasicParticle* copy = new AliFemtoESEBasicParticle(sqrt(pow(particle->P(),2)+pow(0.13957, 2)),particle->Px(),particle->Py(),particle->Pz(),particle->Charge(),particle->Phi(),particle->Eta());
1466 copy->SetPsiEP(psi);
1467 copy->SetTPCClusterMap(particle->GetTPCClusterMap());
1468 copy->SetTPCSharedMap(particle->GetTPCSharedMap());
698fcf10 1469
e7199f47 1470 tracksClone->Add(copy);
698fcf10 1471 }
1472
1473 return tracksClone;
e7199f47 1474 }*/
698fcf10 1475
1476Double_t AliAnalysisTaskFemtoESE::GetDeltaPhiEP(Double_t px1, Double_t py1, Double_t px2, Double_t py2, Double_t psi)
1477{
1478 // angle of pair
1479 Double_t px = px1+px2;
1480 Double_t py = py1+py2;
1481 Double_t phi = atan2(py,px);
1482
1483 Double_t dphi = phi-psi;
d0570384 1484 while(dphi < 0) dphi += 2*TMath::Pi();
1485 while(dphi > 2*TMath::Pi()) dphi -= 2*TMath::Pi();
698fcf10 1486
1487 return dphi;
1488}
1489
aa5e3401 1490Bool_t AliAnalysisTaskFemtoESE::FindBin(Double_t kt, Double_t phi, Double_t cent, Int_t& a, Int_t& b, Int_t&c)
698fcf10 1491{
aa5e3401 1492 a = b = c = -1;
698fcf10 1493 for(Int_t i = 0; i < nKtBins; i++)
1494 {
1495 if(kt >= ktBins[i] && kt < ktBins[i+1])
1496 {
1497 a = i;
1498 break;
1499 }
1500 }
1501 if(a==-1) return kFALSE;
1502
d0570384 1503 if(phi > TMath::Pi()) phi = 2*TMath::Pi()-phi;
698fcf10 1504 for(Int_t i = 0; i < nEPBins; i++)
1505 {
1506 if(phi >= epBins[i] && phi < epBins[i+1])
1507 {
1508 b = i;
1509 break;
1510 }
1511 }
1512 if(b==-1) return kFALSE;
1513
1514 for(Int_t i = 0; i < nCentBins; i++)
1515 {
1516 if(cent >= centBins[i] && cent < centBins[i+1])
1517 {
1518 c = i;
1519 break;
1520 }
1521 }
1522 if(c==-1) return kFALSE;
1523
698fcf10 1524 return kTRUE;
1525}
1526
1527
1528
1529void AliAnalysisTaskFemtoESE::SetKtBins(Int_t n, Double_t* bins)
1530{
1531 if(ktBins) delete [] ktBins;
1532 nKtBins = n;
1533 nKtBins1 = n+1;
1534 ktBins = new Double_t[nKtBins+1];
1535 for(Int_t i = 0; i < nKtBins+1; i++)
1536 ktBins[i]=bins[i];
537c8bdb 1537 Printf("Setting %i kt bins: ",nKtBins);
1538 for(Int_t i = 0; i < nKtBins+1; i++) Printf("%lf",ktBins[i]);
698fcf10 1539}
1540void AliAnalysisTaskFemtoESE::SetCentBins(Int_t n, Double_t* bins)
1541{
1542 if(centBins) delete [] centBins;
1543 nCentBins = n;
1544 nCentBins1 = n+1;
1545 centBins = new Double_t[nCentBins+1];
1546 for(Int_t i = 0; i < nCentBins+1; i++)
1547 centBins[i]=bins[i];
537c8bdb 1548 Printf("Setting %i centrality bins: ",nCentBins);
1549 for(Int_t i = 0; i < nCentBins+1; i++) Printf("%lf",centBins[i]);
698fcf10 1550}
1551void AliAnalysisTaskFemtoESE::SetVzBins(Int_t n, Double_t* bins)
1552{
1553 if(vzBins) delete [] vzBins;
1554 nVzBins = n;
1555 nVzBins1 = n+1;
1556 vzBins = new Double_t[nVzBins+1];
1557 for(Int_t i = 0; i < nVzBins+1; i++)
1558 vzBins[i]=bins[i];
537c8bdb 1559 Printf("Setting %i vz bins: ",nVzBins);
1560 for(Int_t i = 0; i < nVzBins+1; i++) Printf("%lf",vzBins[i]);
698fcf10 1561}
d0570384 1562void AliAnalysisTaskFemtoESE::SetEPBins(Int_t n)
698fcf10 1563{
1564 if(epBins) delete [] epBins;
1565 nEPBins = n;
1566 nEPBins1 = n+1;
1567 epBins = new Double_t[nEPBins+1];
1568 for(Int_t y = 0; y < nEPBins+1; y++)
d0570384 1569 epBins[y] = (TMath::Pi()/(Double_t)n)*((Double_t)y);
537c8bdb 1570 Printf("Setting %i EP bins: ",nEPBins);
1571 for(Int_t i = 0; i < nEPBins+1; i++) Printf("%lf",epBins[i]);
698fcf10 1572}
1573
d0570384 1574/*Double_t AliAnalysisTaskFemtoESE::GetQPercLHC11h(Double_t qvec)
76082c53 1575{
1576 // preliminary attemp at calculating qvector percentile in LHC11h -- still very approximate and only works in 5% bins
1577 if(!qPercBinsLHC11h)
1578 {
1579 Double_t tempArray[21] = {0.0, 0.22995, 0.33047, 0.410831, 0.480728, 0.545566, 0.606841, 0.66634, 0.725193, 0.783813, 0.843311, 0.904185, 0.96796, 1.03522, 1.10768, 1.18774, 1.27808, 1.3857, 1.52438, 1.73633, 4.95};
1580 nqPercBinsLHC11h = 21;
1581 qPercBinsLHC11h = new Double_t[nqPercBinsLHC11h];
1582 for(Int_t n = 0; n < nqPercBinsLHC11h; n++) qPercBinsLHC11h[n] = tempArray[n];
1583 }
1584
1585 for(Int_t t = 0; t < nqPercBinsLHC11h-1; t++)
1586 if(qvec > qPercBinsLHC11h[t] && qvec < qPercBinsLHC11h[t+1]) return 50.*(2*t+1)/(Double_t)(nqPercBinsLHC11h-1);
1587
1588 if(qvec < qPercBinsLHC11h[0]) return 0.0;
1589 if(qvec > qPercBinsLHC11h[nqPercBinsLHC11h-1]) return 100.0;
1590
1591 return 0.0;
1592
d0570384 1593 }*/
f8caeea5 1594
3b2ee55b 1595/*Double_t AliAnalysisTaskFemtoESE::GetCentralityWeight(Double_t cent)
f8caeea5 1596{
1597
1598 // use makeCentWeighting.C to fit centrality distribution to obtain this parameterization
1599 Double_t par1 = 2.60629;
1600 Double_t par2 = 0.579333;
1601 Double_t limit1 = 8.71488;
1602 Double_t limit2 = 12.0126;
1603
1604 if(cent < limit1) return 1./par1;
1605 if(cent > limit2) return 1./par2;
1606
1607 Double_t slope = (par2-par1)/(limit2-limit1);
1608 Double_t b = par1-slope*limit1;
1609
1610 return 1./(slope*cent+b);
3b2ee55b 1611 }*/
1612
1613Double_t AliAnalysisTaskFemtoESE::GetCentralityWeight(Double_t cent)
1614{
1615
1616 // use makeCentWeighting.C to fit centrality distribution to obtain this parameterization
1617 Double_t par1 = 2.60629;
1618 Double_t par2 = 0.579333;
1619 Double_t limit1 = 8.71488;
1620 Double_t limit2 = 12.0126;
1621
1622 if(cent < limit1) return 1.;
1623 if(cent > limit2) return 1.;
1624
1625 Double_t slope = (par2-par1)/(limit2-limit1);
1626 Double_t b = par1-slope*limit1;
1627
1628 if(cent < 10.)
1629 return par1/(slope*cent+b);
1630 else
1631 return par2/(slope*cent+b);
f8caeea5 1632}