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