]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliAnalysisTaskHFEFlow.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add fast merging option (Diego)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskHFEFlow.cxx
CommitLineData
8c1c76e9 1/**************************************************************************\r
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3* *\r
4* Author: The ALICE Off-line Project. *\r
5* Contributors are mentioned in the code where appropriate. *\r
6* *\r
7* Permission to use, copy, modify and distribute this software and its *\r
8* documentation strictly for non-commercial purposes is hereby granted *\r
9* without fee, provided that the above copyright notice appears in all *\r
10* copies and that both the copyright notice and this permission notice *\r
11* appear in the supporting documentation. The authors make no claims *\r
12* about the suitability of this software for any purpose. It is *\r
13* provided "as is" without express or implied warranty. *\r
14**************************************************************************/\r
15//\r
16// Flow task\r
17// \r
18// Authors:\r
19// Raphaelle Bailhache <R.Bailhache@gsi.de>\r
20//\r
21#include "TROOT.h"\r
22#include "TH1D.h"\r
23#include "TH2D.h"\r
24#include "TChain.h"\r
25#include "TVector2.h"\r
26#include "THnSparse.h"\r
27#include "TMath.h"\r
28#include "TRandom3.h"\r
29\r
30#include "AliAnalysisTaskSE.h"\r
31\r
32#include "AliESDInputHandler.h"\r
33#include "AliMCEvent.h"\r
34#include "AliESD.h"\r
35#include "AliESDEvent.h"\r
36#include "AliPIDResponse.h"\r
37#include "AliESDVZERO.h"\r
38#include "AliESDUtils.h"\r
39#include "AliMCParticle.h"\r
40\r
41#include "AliFlowCandidateTrack.h"\r
42#include "AliFlowEvent.h"\r
43#include "AliFlowTrackCuts.h"\r
44#include "AliFlowVector.h"\r
45#include "AliFlowCommonConstants.h"\r
46\r
47\r
48#include "AliHFEcuts.h"\r
49#include "AliHFEpid.h"\r
50#include "AliHFEpidQAmanager.h"\r
51#include "AliHFEtools.h"\r
52\r
53\r
54\r
55#include "AliCentrality.h"\r
56#include "AliEventplane.h"\r
57#include "AliAnalysisTaskHFEFlow.h"\r
58\r
59\r
60//____________________________________________________________________\r
61AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :\r
62 AliAnalysisTaskSE(),\r
63 fListHist(), \r
64 fSubEtaGapTPC(kFALSE),\r
65 fEtaGap(0.0),\r
66 fNbBinsPtQCumulant(15),\r
67 fMinPtQCumulant(0.0),\r
68 fMaxPtQCumulant(6.0),\r
69 fAfterBurnerOn(kFALSE),\r
70 fNonFlowNumberOfTrackClones(0),\r
71 fV1(0.),\r
72 fV2(0.),\r
73 fV3(0.),\r
74 fV4(0.),\r
75 fV5(0.),\r
76 fMaxNumberOfIterations(100),\r
77 fPrecisionPhi(0.001),\r
78 fUseMCReactionPlane(kFALSE),\r
79 fMCPID(kFALSE),\r
80 fcutsRP(0),\r
81 fcutsPOI(0),\r
82 fHFECuts(0),\r
83 fPID(0),\r
84 fPIDqa(0),\r
85 fflowEvent(0x0),\r
86 fHistEV(0),\r
87 fEventPlane(0x0),\r
88 fEventPlaneaftersubtraction(0x0),\r
89 fCos2phie(0x0),\r
90 fSin2phie(0x0),\r
91 fCos2phiep(0x0),\r
92 fSin2phiep(0x0),\r
93 fSin2phiephiep(0x0),\r
94 fCosRes(0x0),\r
95 fDeltaPhiMaps(0x0),\r
96 fCosPhiMaps(0x0)\r
97{\r
98 // Constructor\r
99\r
100 \r
101 \r
102}\r
103//______________________________________________________________________________\r
104AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :\r
105 AliAnalysisTaskSE(name),\r
106 fListHist(), \r
107 fSubEtaGapTPC(kFALSE),\r
108 fEtaGap(0.0),\r
109 fNbBinsPtQCumulant(15),\r
110 fMinPtQCumulant(0.0),\r
111 fMaxPtQCumulant(6.0),\r
112 fAfterBurnerOn(kFALSE),\r
113 fNonFlowNumberOfTrackClones(0),\r
114 fV1(0.),\r
115 fV2(0.),\r
116 fV3(0.),\r
117 fV4(0.),\r
118 fV5(0.),\r
119 fMaxNumberOfIterations(100),\r
120 fPrecisionPhi(0.001),\r
121 fUseMCReactionPlane(kFALSE),\r
122 fMCPID(kFALSE),\r
123 fcutsRP(0),\r
124 fcutsPOI(0),\r
125 fHFECuts(0),\r
126 fPID(0),\r
127 fPIDqa(0),\r
128 fflowEvent(0x0),\r
129 fHistEV(0),\r
130 fEventPlane(0x0),\r
131 fEventPlaneaftersubtraction(0x0),\r
132 fCos2phie(0x0),\r
133 fSin2phie(0x0),\r
134 fCos2phiep(0x0),\r
135 fSin2phiep(0x0),\r
136 fSin2phiephiep(0x0),\r
137 fCosRes(0x0),\r
138 fDeltaPhiMaps(0x0),\r
139 fCosPhiMaps(0x0)\r
140{\r
141 //\r
142 // named ctor\r
143 //\r
144\r
145 fPID = new AliHFEpid("hfePid");\r
146 fPIDqa = new AliHFEpidQAmanager;\r
147\r
148 DefineInput(0,TChain::Class());\r
149 DefineOutput(1, TList::Class());\r
150 DefineOutput(2,AliFlowEventSimple::Class()); // first band 0-20%\r
151 DefineOutput(3,AliFlowEventSimple::Class()); // first band 20-40%\r
152 DefineOutput(4,AliFlowEventSimple::Class()); // first band 40-60%\r
153 DefineOutput(5,AliFlowEventSimple::Class()); // first band 60-80%\r
154 \r
155}\r
156//________________________________________________________________________\r
157void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()\r
158{\r
159\r
160 //********************\r
161 // Create histograms\r
162 //********************\r
163\r
164 //**************\r
165 // Cuts\r
166 //**************\r
167\r
168 //---------Data selection----------\r
169 //kMC, kGlobal, kTPCstandalone, kSPDtracklet, kPMD\r
170 //AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;\r
171 //AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;\r
172\r
173 //---------Parameter mixing--------\r
174 //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt\r
175 //AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;\r
176 //AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;\r
177\r
178 // RP TRACK CUTS:\r
179 fcutsRP = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();\r
180 fcutsRP->SetName("StandartTPC");\r
181 fcutsRP->SetEtaRange(-0.9,0.9);\r
182 fcutsRP->SetQA(kTRUE);\r
183\r
184 //POI TRACK CUTS:\r
185 fcutsPOI = new AliFlowTrackCuts("dummy");\r
186 fcutsPOI->SetParamType(AliFlowTrackCuts::kGlobal);\r
187 fcutsPOI->SetPtRange(+1,-1); // select nothing QUICK\r
188 fcutsPOI->SetEtaRange(+1,-1); // select nothing VZERO\r
189\r
190 // Flow\r
191 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();\r
192 cc->SetNbinsMult(10000);\r
193 cc->SetMultMin(0);\r
194 cc->SetMultMax(10000.);\r
195 cc->SetNbinsPt(fNbBinsPtQCumulant);\r
196 cc->SetPtMin(fMinPtQCumulant);\r
197 cc->SetPtMax(fMaxPtQCumulant);\r
198 cc->SetNbinsPhi(180);\r
199 cc->SetPhiMin(0.0);\r
200 cc->SetPhiMax(TMath::TwoPi());\r
201 cc->SetNbinsEta(200);\r
202 cc->SetEtaMin(-0.9);\r
203 cc->SetEtaMax(+0.9);\r
204 cc->SetNbinsQ(500);\r
205 cc->SetQMin(0.0);\r
206 cc->SetQMax(3.0);\r
207\r
208 \r
209 // HFE cuts\r
210\r
211 if(!fHFECuts){\r
212 fHFECuts = new AliHFEcuts;\r
213 fHFECuts->CreateStandardCuts();\r
214 }\r
215 fHFECuts->Initialize();\r
216 \r
217 // PID HFE\r
218 //fPID->SetHasMCData(HasMCData());\r
219 if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);\r
220 fPID->InitializePID();\r
221 fPIDqa->Initialize(fPID);\r
222 fPID->SortDetectors();\r
223 \r
224 //**************************\r
225 // Bins for the THnSparse\r
226 //**************************\r
227\r
228 Int_t nBinsPt = 25;\r
229 Double_t minPt = 0.001;\r
230 Double_t maxPt = 10.0;\r
231 Double_t binLimLogPt[nBinsPt+1];\r
232 Double_t binLimPt[nBinsPt+1];\r
233 for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;\r
234 for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);\r
235\r
236 Int_t nBinsEta = 8;\r
237 Double_t minEta = -0.8;\r
238 Double_t maxEta = 0.8;\r
239 Double_t binLimEta[nBinsEta+1];\r
240 for(Int_t i=0; i<=nBinsEta; i++) binLimEta[i]=(Double_t)minEta + (maxEta-minEta)/nBinsEta*(Double_t)i ;\r
241 \r
242 Int_t nBinsCos = 50;\r
243 Double_t minCos = -1.0;\r
244 Double_t maxCos = 1.0;\r
245 Double_t binLimCos[nBinsCos+1];\r
246 for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;\r
247 \r
248\r
249 Int_t nBinsC = 11;\r
250 Double_t minC = 0.0;\r
251 Double_t maxC = 11.0;\r
252 Double_t binLimC[nBinsC+1];\r
253 for(Int_t i=0; i<=nBinsC; i++) binLimC[i]=(Double_t)minC + (maxC-minC)/nBinsC*(Double_t)i ;\r
254 \r
255 Int_t nBinsPhi = 50;\r
256 Double_t minPhi = 0.0;\r
257 Double_t maxPhi = TMath::Pi();\r
258 Double_t binLimPhi[nBinsPhi+1];\r
259 for(Int_t i=0; i<=nBinsPhi; i++) {\r
260 binLimPhi[i]=(Double_t)minPhi + (maxPhi-minPhi)/nBinsPhi*(Double_t)i ;\r
261 //printf("bin phi is %f for %d\n",binLimPhi[i],i);\r
262 }\r
263 \r
264 Int_t nBinsQ = 200;\r
265 Double_t minQ = -300.0;\r
266 Double_t maxQ = 300.0;\r
267 Double_t binLimQ[nBinsQ+1];\r
268 for(Int_t i=0; i<=nBinsQ; i++) binLimQ[i]=(Double_t)minQ + (maxQ-minQ)/nBinsQ*(Double_t)i ;\r
269\r
270 Int_t nBinsM = 250;\r
271 Double_t minM = 0.0;\r
272 Double_t maxM = 800.0;\r
273 Double_t binLimM[nBinsM+1];\r
274 for(Int_t i=0; i<=nBinsM; i++) binLimM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;\r
275\r
276 Int_t nBinsCELL = 64;\r
277 Double_t minCELL = 0.0;\r
278 Double_t maxCELL = 64.0;\r
279 Double_t binLimCELL[nBinsCELL+1];\r
280 for(Int_t i=0; i<=nBinsCELL; i++) binLimCELL[i]=(Double_t)minCELL + (maxCELL-minCELL)/nBinsCELL*(Double_t)i ;\r
281 \r
282\r
283 //******************\r
284 // Histograms\r
285 //******************\r
286 \r
287 fListHist = new TList();\r
288\r
289 // Histos\r
290 fHistEV = new TH2D("fHistEV", "events", 3, 0, 3, 3, 0,3);\r
291 \r
292 // Event plane as function of phiep, centrality\r
293 const Int_t nDima=2;\r
294 Int_t nBina[nDima] = {nBinsPhi, nBinsC};\r
295 fEventPlane = new THnSparseF("EventPlane","EventPlane",nDima,nBina);\r
296 fEventPlane->SetBinEdges(0,binLimPhi);\r
297 fEventPlane->SetBinEdges(1,binLimC);\r
298 \r
299 // Event Plane after subtraction as function of phiep, centrality, pt, eta\r
300 const Int_t nDimb=2;\r
301 Int_t nBinb[nDimb] = {nBinsPhi, nBinsC};\r
302 fEventPlaneaftersubtraction = new THnSparseF("EventPlane_aftersubtraction","EventPlane_aftersubtraction",nDimb,nBinb);\r
303 fEventPlaneaftersubtraction->SetBinEdges(0,binLimPhi);\r
304 fEventPlaneaftersubtraction->SetBinEdges(1,binLimC);\r
305\r
306 // Monitoring Event plane after subtraction of the track\r
307 const Int_t nDime=4;\r
308 Int_t nBine[nDime] = {nBinsCos, nBinsC, nBinsPt, nBinsEta};\r
309 fCos2phie = new THnSparseF("cos2phie","cos2phie",nDime,nBine);\r
310 fCos2phie->SetBinEdges(2,binLimPt);\r
311 fCos2phie->SetBinEdges(3,binLimEta);\r
312 fCos2phie->SetBinEdges(0,binLimCos);\r
313 fCos2phie->SetBinEdges(1,binLimC);\r
314 fSin2phie = new THnSparseF("sin2phie","sin2phie",nDime,nBine);\r
315 fSin2phie->SetBinEdges(2,binLimPt);\r
316 fSin2phie->SetBinEdges(3,binLimEta);\r
317 fSin2phie->SetBinEdges(0,binLimCos);\r
318 fSin2phie->SetBinEdges(1,binLimC);\r
319 fCos2phiep = new THnSparseF("cos2phiep","cos2phiep",nDime,nBine);\r
320 fCos2phiep->SetBinEdges(2,binLimPt);\r
321 fCos2phiep->SetBinEdges(3,binLimEta);\r
322 fCos2phiep->SetBinEdges(0,binLimCos);\r
323 fCos2phiep->SetBinEdges(1,binLimC);\r
324 fSin2phiep = new THnSparseF("sin2phiep","sin2phiep",nDime,nBine);\r
325 fSin2phiep->SetBinEdges(2,binLimPt);\r
326 fSin2phiep->SetBinEdges(3,binLimEta);\r
327 fSin2phiep->SetBinEdges(0,binLimCos);\r
328 fSin2phiep->SetBinEdges(1,binLimC);\r
329 fSin2phiephiep = new THnSparseF("sin2phie_phiep","sin2phie_phiep",nDime,nBine);\r
330 fSin2phiephiep->SetBinEdges(2,binLimPt);\r
331 fSin2phiephiep->SetBinEdges(3,binLimEta);\r
332 fSin2phiephiep->SetBinEdges(0,binLimCos);\r
333 fSin2phiephiep->SetBinEdges(1,binLimC);\r
334 \r
335 // Resolution cosres centrality\r
336 const Int_t nDimf=2;\r
337 Int_t nBinf[nDimf] = {nBinsCos, nBinsC};\r
338 fCosRes = new THnSparseF("CosRes","CosRes",nDimf,nBinf);\r
339 fCosRes->SetBinEdges(0,binLimCos);\r
340 fCosRes->SetBinEdges(1,binLimC);\r
341 \r
342 // Maps delta phi\r
343 const Int_t nDimg=3;\r
344 Int_t nBing[nDimg] = {nBinsPhi,nBinsC,nBinsPt};\r
345 fDeltaPhiMaps = new THnSparseF("DeltaPhiMaps","DeltaPhiMaps",nDimg,nBing);\r
346 fDeltaPhiMaps->SetBinEdges(0,binLimPhi);\r
347 fDeltaPhiMaps->SetBinEdges(1,binLimC);\r
348 fDeltaPhiMaps->SetBinEdges(2,binLimPt);\r
349 \r
350 // Maps cos phi\r
351 const Int_t nDimh=3;\r
352 Int_t nBinh[nDimh] = {nBinsCos,nBinsC,nBinsPt};\r
353 fCosPhiMaps = new THnSparseF("CosPhiMaps","CosPhiMaps",nDimh,nBinh);\r
354 fCosPhiMaps->SetBinEdges(0,binLimCos);\r
355 fCosPhiMaps->SetBinEdges(1,binLimC);\r
356 fCosPhiMaps->SetBinEdges(2,binLimPt);\r
357 \r
358 //**************************\r
359 // Add to the list\r
360 //******************************\r
361\r
362 fListHist->Add(fPIDqa->MakeList("HFEpidQA"));\r
363 \r
364 fListHist->Add(fHistEV);\r
365\r
366 fListHist->Add(fEventPlane);\r
367 fListHist->Add(fEventPlaneaftersubtraction);\r
368 fListHist->Add(fCos2phie);\r
369 fListHist->Add(fSin2phie);\r
370 fListHist->Add(fCos2phiep);\r
371 fListHist->Add(fSin2phiep);\r
372 fListHist->Add(fSin2phiephiep);\r
373 fListHist->Add(fCosRes);\r
374 fListHist->Add(fDeltaPhiMaps);\r
375 fListHist->Add(fCosPhiMaps);\r
376 \r
377\r
378 PostData(1, fListHist);\r
379\r
380\r
381}\r
382 \r
383//________________________________________________________________________\r
384void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)\r
385{\r
386 //\r
387 // Loop over event\r
388 //\r
389 \r
390 Double_t massElectron = 0.000511;\r
391 Double_t mcReactionPlane = 0.0;\r
392\r
393 Float_t cntr = 0.0;\r
394 Double_t binct = 11.5;\r
395 Float_t binctt = -1.0;\r
396 \r
397 Double_t valuensparsea[2];\r
398 Double_t valuensparsee[4];\r
399 Double_t valuensparsef[2];\r
400 Double_t valuensparseg[3];\r
401 Double_t valuensparseh[3];\r
402\r
403 AliMCEvent *mcEvent = MCEvent();\r
404 AliMCParticle *mctrack = NULL;\r
405 \r
406 /////////////////\r
407 // centrality\r
408 /////////////////\r
409 \r
410 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
411 if(!esd) return;\r
412 AliCentrality *centrality = esd->GetCentrality();\r
413 if(!centrality) return;\r
414 cntr = centrality->GetCentralityPercentile("V0M");\r
415 if((0.0< cntr) && (cntr<5.0)) binct = 0.5;\r
416 if((5.0< cntr) && (cntr<10.0)) binct = 1.5;\r
417 if((10.0< cntr) && (cntr<20.0)) binct = 2.5;\r
418 if((20.0< cntr) && (cntr<30.0)) binct = 3.5;\r
419 if((30.0< cntr) && (cntr<40.0)) binct = 4.5;\r
420 if((40.0< cntr) && (cntr<50.0)) binct = 5.5;\r
421 if((50.0< cntr) && (cntr<60.0)) binct = 6.5;\r
422 if((60.0< cntr) && (cntr<70.0)) binct = 7.5;\r
423 if((70.0< cntr) && (cntr<80.0)) binct = 8.5;\r
424 if((80.0< cntr) && (cntr<90.0)) binct = 9.5;\r
425 if((90.0< cntr) && (cntr<100.0)) binct = 10.5;\r
426 \r
427 if((0.< cntr) && (cntr < 20.)) binctt = 0.5;\r
428 if((20.< cntr) && (cntr < 40.)) binctt = 1.5;\r
429 if((40.< cntr) && (cntr < 80.)) binctt = 2.5;\r
430\r
431 if(binct > 11.0) return;\r
432 \r
433 // centrality\r
434 valuensparsea[1] = binct; \r
435 valuensparsee[1] = binct; \r
436 valuensparsef[1] = binct; \r
437 valuensparseg[1] = binct;\r
438 valuensparseh[1] = binct; \r
439 \r
440\r
441 //////////////////////\r
442 // run number\r
443 //////////////////////\r
444\r
445 Int_t runnumber = esd->GetRunNumber();\r
446 \r
447 if(!fPID->IsInitialized()){\r
448 // Initialize PID with the given run number\r
449 fPID->InitializePID(runnumber);\r
450 }\r
451\r
452\r
453 //////////\r
454 // PID\r
455 //////////\r
456 \r
457 AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();\r
458 if(!pidResponse){\r
459 //printf("No PID response set\n");\r
460 return;\r
461 }\r
462 fPID->SetPIDResponse(pidResponse);\r
463\r
464 fHistEV->Fill(binctt,0.0);\r
465 \r
466\r
467 //////////////////\r
468 // Event cut\r
469 //////////////////\r
470 if(!fHFECuts->CheckEventCuts("fEvRecCuts", esd)) {\r
471 PostData(1, fListHist);\r
472 return;\r
473 }\r
474\r
475 fHistEV->Fill(binctt,1.0);\r
476\r
477 //////////////////////////////////// \r
478 // First method TPC event plane\r
479 ////////////////////////////////////\r
480\r
481 AliEventplane* esdEPa = esd->GetEventplane();\r
482 TVector2 *standardQ = esdEPa->GetQVector(); // This is the "standard" Q-Vector\r
483 Double_t qx = -1.0;\r
484 Double_t qy = -1.0;\r
485 if(standardQ) {\r
486 qx = standardQ->X();\r
487 qy = standardQ->Y();\r
488 } \r
489 TVector2 qVectorfortrack;\r
490 qVectorfortrack.Set(qx,qy);\r
491 Float_t eventPlanea = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.; \r
492 \r
493 TVector2 *qsub1a = esdEPa->GetQsub1();\r
494 TVector2 *qsub2a = esdEPa->GetQsub2();\r
495 Float_t eventPlanesub1a = -100.0;\r
496 Float_t eventPlanesub2a = -100.0;\r
497 if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;\r
498 if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;\r
499 Double_t diffsub1sub2a = -100.0;\r
500 if(qsub1a && qsub2a) {\r
501 diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
502 }\r
503\r
504 //////////////////////////////////////////\r
505 // Cut for event with TPC event defined\r
506 /////////////////////////////////////////\r
507 \r
508 if(!standardQ) {\r
509 PostData(1, fListHist);\r
510 return;\r
511 }\r
512\r
513 ///////////////////////\r
514 // AliFlowEvent \r
515 //////////////////////\r
516\r
517 Int_t nbtracks = esd->GetNumberOfTracks();\r
518 //printf("Number of tracks %d\n",nbtracks);\r
519\r
520 fcutsRP->SetEvent( InputEvent(), MCEvent());\r
521 fcutsPOI->SetEvent( InputEvent(), MCEvent());\r
522 if( fflowEvent ){ \r
523 fflowEvent->~AliFlowEvent();\r
524 new(fflowEvent) AliFlowEvent(fcutsRP,fcutsPOI);\r
525 }else fflowEvent = new AliFlowEvent(fcutsRP,fcutsPOI);\r
526 if(mcEvent && mcEvent->GenEventHeader()) {\r
527 fflowEvent->SetMCReactionPlaneAngle(mcEvent);\r
528 //if reaction plane not set from elsewhere randomize it before adding flow\r
529 //if (!fflowEvent->IsSetMCReactionPlaneAngle()) fflowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));\r
530 mcReactionPlane = TVector2::Phi_0_2pi(fflowEvent->GetMCReactionPlaneAngle());\r
531 if(mcReactionPlane > TMath::Pi()) mcReactionPlane = mcReactionPlane - TMath::Pi();\r
532 //printf("MC reaction plane %f\n",mcReactionPlane);\r
533 }\r
534 fflowEvent->SetReferenceMultiplicity( nbtracks );\r
535 fflowEvent->DefineDeadZone(0,0,0,0);\r
536 //fflowEvent.TagSubeventsInEta(-0.8,-0.1,0.1,0.8);\r
537\r
538 ////////////////\r
539 // MC\r
540 ///////////////\r
541 if(fUseMCReactionPlane) {\r
542 eventPlanea = mcReactionPlane;\r
543 diffsub1sub2a = 0.0;\r
544 }\r
545\r
546 \r
547 //////////////////////\r
548 // Fill Histos\r
549 //////////////////////\r
550\r
551 fHistEV->Fill(binctt,2.0);\r
552\r
553 // Fill\r
554 valuensparsea[0] = eventPlanea; \r
555 fEventPlane->Fill(&valuensparsea[0]);\r
556\r
557 valuensparsef[0] = diffsub1sub2a;\r
558 fCosRes->Fill(&valuensparsef[0]);\r
559 \r
560 //////////////////////////\r
561 // Loop over ESD track\r
562 //////////////////////////\r
563 \r
564\r
565 for(Int_t k = 0; k < nbtracks; k++){\r
566 \r
567 AliESDtrack *track = esd->GetTrack(k);\r
568 if(!track) continue;\r
569\r
570 Bool_t survived = kTRUE;\r
571 for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut <= AliHFEcuts::kStepHFEcutsTRD; icut++){\r
572 if(!fHFECuts->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){\r
573 survived = kFALSE;\r
574 break;\r
575 }\r
576 //printf("Pass the cut %d\n",icut);\r
577 }\r
578 \r
579 if(!survived) continue;\r
580\r
581 // Apply PID for Data\r
582 if(!fMCPID) {\r
583 AliHFEpidObject hfetrack;\r
584 hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
585 hfetrack.SetRecTrack(track);\r
586 hfetrack.SetCentrality((Int_t)binct);\r
587 //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());\r
588 hfetrack.SetPbPb();\r
589 if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) {\r
590 continue;\r
591 }\r
592 }\r
593 else {\r
594 if(!mcEvent) continue;\r
595 if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;\r
596 //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));\r
597 if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;\r
598 }\r
599\r
600\r
601 /////////////////////////////////////////////////////////////////////////////\r
602 // Add candidate to AliFlowEvent for POI and subtract from RP if needed\r
603 ////////////////////////////////////////////////////////////////////////////\r
604 Int_t idtrack = static_cast<AliVTrack*>(track)->GetID();\r
605 Bool_t found = kFALSE;\r
606 Int_t numberoffound = 0;\r
607 //printf("A: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
608 for(Int_t iRPs=0; iRPs< fflowEvent->NumberOfTracks(); iRPs++) {\r
609 AliFlowTrack *iRP = (AliFlowTrack*) (fflowEvent->GetTrack(iRPs));\r
610 //if(!iRP->InRPSelection()) continue;\r
611 if( TMath::Abs(idtrack) == TMath::Abs(iRP->GetID()) ) {\r
612 iRP->SetForPOISelection(kTRUE);\r
613 found = kTRUE;\r
614 numberoffound ++;\r
615 }\r
616 }\r
617 //printf("Found %d mal\n",numberoffound);\r
618 if(!found) {\r
619 AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*) MakeTrack(massElectron,track->Pt(),track->Phi(), track->Eta());\r
620 sTrack->SetID(idtrack);\r
621 fflowEvent->AddTrack(sTrack);\r
622 //printf("Add the track\n");\r
623 }\r
624 //printf("B: Number of tracks %d\n",fflowEvent->NumberOfTracks());\r
625 \r
626 \r
627 /////////////////////////////////////////////////////////\r
628 // Subtract electron candidate from TPC event plane\r
629 ////////////////////////////////////////////////////////\r
630\r
631 // Subtract the tracks from the event plane\r
632 Double_t qX = standardQ->X() - esdEPa->GetQContributionX(track); //Modify the components: subtract the track you want to look at with your analysis\r
633 Double_t qY = standardQ->Y() - esdEPa->GetQContributionY(track); //Modify the components: subtract the track you want to look at with your analysis\r
634 TVector2 newQVectorfortrack;\r
635 newQVectorfortrack.Set(qX,qY);\r
636 Float_t eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; \r
637\r
638 ////////////////////////////////////////\r
639 // Fill pt and eta for the THnSparseF\r
640 ///////////////////////////////////////\r
641\r
642 valuensparsee[2] = track->Pt();\r
643 valuensparsee[3] = track->Eta(); \r
644 valuensparseg[2] = track->Pt();\r
645 valuensparseh[2] = track->Pt();\r
646\r
647 Bool_t fillTPC = kTRUE;\r
648 if((!qsub1a) || (!qsub2a)) fillTPC = kFALSE;\r
649\r
650 if(fSubEtaGapTPC) {\r
651 if(track->Eta() < (- fEtaGap/2.)) eventplanesubtracted = eventPlanesub1a;\r
652 else if(track->Eta() > (fEtaGap/2.)) eventplanesubtracted = eventPlanesub2a;\r
653 else fillTPC = kFALSE;\r
654 }\r
655\r
656 ///////////////\r
657 // MC\r
658 //////////////\r
659 if(fUseMCReactionPlane) {\r
660 eventplanesubtracted = mcReactionPlane;\r
661 fillTPC = kTRUE;\r
662 }\r
663 \r
664 //////////////////////////////////////////////////////////////////////////////\r
665 ///////////////////////////AFTERBURNER\r
666 Double_t phitrack = track->Phi(); \r
667 if (fAfterBurnerOn)\r
668 {\r
669 phitrack = GetPhiAfterAddV2(track->Phi(),mcReactionPlane);\r
670 }\r
671 //////////////////////////////////////////////////////////////////////////////\r
672\r
673\r
674 ///////////////////////\r
675 // Calculate deltaphi\r
676 ///////////////////////\r
677 \r
678 // Suppose phi track is between 0.0 and phi\r
679 Double_t deltaphi = TVector2::Phi_0_2pi(phitrack - eventplanesubtracted);\r
680 if(deltaphi > TMath::Pi()) deltaphi = deltaphi - TMath::Pi();\r
681 \r
682 /////////////////////\r
683 // Fill THnSparseF\r
684 /////////////////////\r
685\r
686 //\r
687 valuensparsea[0] = eventplanesubtracted;\r
688 if(fillTPC) fEventPlaneaftersubtraction->Fill(&valuensparsea[0]);\r
689 \r
690 //\r
691 valuensparsee[0] = TMath::Cos(2*phitrack);\r
692 fCos2phie->Fill(&valuensparsee[0]);\r
693 valuensparsee[0] = TMath::Sin(2*phitrack);\r
694 fSin2phie->Fill(&valuensparsee[0]);\r
695 \r
696 //\r
697 valuensparsee[0] = TMath::Cos(2*eventplanesubtracted);\r
698 if(fillTPC) fCos2phiep->Fill(&valuensparsee[0]);\r
699 valuensparsee[0] = TMath::Sin(2*eventplanesubtracted);\r
700 if(fillTPC) fSin2phiep->Fill(&valuensparsee[0]);\r
701 valuensparsee[0] = TMath::Sin(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
702 if(fillTPC) fSin2phiephiep->Fill(&valuensparsee[0]);\r
703 \r
704 // \r
705 valuensparseg[0] = deltaphi;\r
706 if(fillTPC) fDeltaPhiMaps->Fill(&valuensparseg[0]);\r
707 \r
708 //\r
709 valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
710 if(fillTPC) fCosPhiMaps->Fill(&valuensparseh[0]);\r
711 \r
712 }\r
713\r
714 //////////////////////////////////////////////////////////////////////////////\r
715 ///////////////////////////AFTERBURNER\r
716 if (fAfterBurnerOn)\r
717 {\r
718 fflowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow\r
719 fflowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks\r
720 }\r
721 //////////////////////////////////////////////////////////////////////////////\r
722\r
723\r
724\r
725 if((0.0< cntr) && (cntr<20.0)) PostData(2,fflowEvent);\r
726 if((20.0< cntr) && (cntr<40.0)) PostData(3,fflowEvent);\r
727 if((40.0< cntr) && (cntr<60.0)) PostData(4,fflowEvent);\r
728 if((60.0< cntr) && (cntr<80.0)) PostData(5,fflowEvent);\r
729 \r
730 \r
731 PostData(1, fListHist);\r
732\r
733\r
734 \r
735}\r
736//______________________________________________________________________________\r
737AliFlowCandidateTrack *AliAnalysisTaskHFEFlow::MakeTrack( Double_t mass, \r
738 Double_t pt, Double_t phi, Double_t eta) {\r
739 //\r
740 // Make Track (Not needed actually)\r
741 //\r
742\r
743 AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();\r
744 sTrack->SetMass(mass);\r
745 sTrack->SetPt(pt);\r
746 sTrack->SetPhi(phi);\r
747 sTrack->SetEta(eta);\r
748 sTrack->SetForPOISelection(kTRUE);\r
749 sTrack->SetForRPSelection(kFALSE);\r
750 return sTrack;\r
751}\r
752//_________________________________________________________________________________ \r
753Double_t AliAnalysisTaskHFEFlow::GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const\r
754{\r
755 //\r
756 // Adds v2, uses Newton-Raphson iteration\r
757 //\r
758 Double_t phiend=phi;\r
759 Double_t phi0=phi;\r
760 Double_t f=0.;\r
761 Double_t fp=0.;\r
762 Double_t phiprev=0.;\r
763\r
764 for (Int_t i=0; i<fMaxNumberOfIterations; i++)\r
765 {\r
766 phiprev=phiend; //store last value for comparison\r
767 f = phiend-phi0+fV2*TMath::Sin(2.*(phiend-reactionPlaneAngle));\r
768 fp = 1.0+2.0*fV2*TMath::Cos(2.*(phiend-reactionPlaneAngle)); //first derivative\r
769 phiend -= f/fp;\r
770 if (TMath::AreEqualAbs(phiprev,phiend,fPrecisionPhi)) break;\r
771 }\r
772 return phiend;\r
773}\r